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ABSTRACT 


Emitter location finding enables important functionality for both military and 
civilian applications. GPS is the most recognized and widely used positioning 
system, but it is a receiver location system that functions in a markedly different 
manner from emitter location. Many geo-location techniques predate and have 
been proposed as an alternative to GPS. Some of the more commonly used and 
exploited of these techniques are angle of arrival, time of arrival, frequency 
difference of arrival, and time difference of arrival (TDOA). This thesis is primarily 
focused on TDOA. 

These techniques are important for military applications. Location finding 
is a part of electronic warfare support, which is one of the main braches of 
electronic warfare. Because these techniques are platform independent, they can 
be used with any system or platform, such as UAVs, manned aircraft, ground 
locations, etc. In Turkey it is vitally important for the army conducting search and 
destroy operations against terrorists to locate emitters associated with these 
terrorists. 

The simulation developed in this thesis provides a better understanding of 
the accuracy of TDOA based geolocation systems. Combinations of receivers 
and techniques are explored to determine the optimal solutions. The factors of 
noise and distance have a linear effect on accuracy. The best combination of 
receivers is determined with consideration to using a combination of fixed and 
airborne platforms. The best distribution for highest accuracy is determined and 
discussed. 
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I. INTRODUCTION 


A. AREA OF RESEARCH 

Electronic warfare (EW) plays a dominant role in today’s world of 
technological warfare. It has been considered to be a force multiplier for 
decades. Nations who understand the importance of the EW have always used it 
and benefited. We as the soldiers cannot imagine a battlefield without EW and its 
components. “EW is one of the five core capabilities” (JP 3-13-1, 2007). Joint 
Publication (JP) 3-13-1 Electronic Warfare is U.S. joint doctrine for EW and 
provides joint doctrine for every part of EW, ranging from planning and 
preparation to execution for military operations. This publication is an important 
source for anyone wanting to define what EW is and what components compose 
it. We cannot win a war or even a battle without the use of EW. 

EW includes three major branches: Electronic Attack (EA), Electronic 
Protection (EP) and Electronic Warfare Support (ES). “EA involves the use of EM 
energy, directed energy, or anti-radiation weapons to attack personnel, facilities, 
or equipment with the intent of degrading, neutralizing, or destroying enemy 
combat capability and is considered a form of fire. EP involves actions taken to 
protect personnel, facilities, and equipment from any effects of friendly or enemy 
use of the electromagnetic spectrum that degrade, neutralize, or destroy friendly 
combat capability. ES is the subdivision of EW involving actions tasked by, or 
under direct control of, an operational commander to search for, intercept, 
identify, and locate or localize sources of intentional and unintentional radiated 
EM energy for the purpose of immediate threat recognition, targeting, planning, 
and conduct of future operations” (JP 3-13-1,2007). 

ES includes direction finding of enemy RF signals. Direction finding can 
also be considered as finding the enemy location by analyzing the signal or 
signals that are transmitted by the emitter. This is also called emitter location or, 
if the emitter is on the earth’s surface, emitter geo-location. 
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It is important for both the military and civilians to be able to find the 
location of RF emitters on the surface of the earth under extreme conditions, 
such as in war for the military and in the case of accidents or disasters for civilian 
response. Military forces want to know the where the enemy is. Determining the 
location of enemy forces has always been important to combat, and various 
methods have been studied and developed through the ages. The easiest and 
most accurate way to do this historically has been to send a man forward to 
locate the enemy. Time and available human resources limit this. There are 
inherent risks, such as having the soldiers killed and being recognized by the 
enemy forces if they capture the reconnaissance team. Because of these 
limitations people have tried to use technology as much as possible in the 
modern era. 

It is important to know the position of the enemy without letting them know 
that you are trying to locate them. This can be done with passive location finding 
techniques. In contrast, people have been using active location finding 
techniques, such as radar, for decades. When active location finding techniques 
are used they let the enemy know that they are being identified or located 
because active location finding techniques uses detectable RF signals. 

To reduce the chance that the enemy will discover that they are under 
scrutiny, scientists and engineers have developed techniques for passive 
location finding. These techniques depend upon the interception of enemy 
transmissions within the EM spectrum, rather than emanating friendly EM 
transmissions. Some of the more common of these techniques are 

• TDOA (Time Difference of Arrival) 

• FDOA (Frequency Difference of Arrival) 

• AOA (Angle of Arrival) 

• TOA (Time of Arrival) 
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Among these techniques, TDOA is currently the most commonly studied 
technique. This technique is platform independent. TDOA is most easily applied 
to pulsed signals. Since all signals of interest are not pulsed signals, some 
techniques have to be implemented for continuous wave (CW) signals. Pulsed 
signals have distinctive features that can be recognized in the waveform. CW 
signals do not have the same distinctive features and can change over time. The 
signal depends on the message that is carried over the carrier. Since the 
message over the carrier is random, it is difficult to determine the beginning and 
the end of the signal. 

Radars typically use pulsed signals, and communication systems typically 
use various modulations of CW signals. Ground forces are normally more 
interested in communication signals than radar; military forces need to coordinate 
their needs with scientists to find an appropriate geolocation technique for the 
CW signals associated with these communications. 

The FDOA technique requires a long-duration signal to be able to 
determine direction to a sufficient accuracy for position determination. 
Conversely, the TDOA technique does not need this long duration. Instead of 
requiring a long duration signal, very short duration signals can be used to 
determine the location if all the necessary receivers receive the signal. Proper 
coordination of the signals between the receivers and the main node that 
compute the position, and synchronization among the receivers, are necessary 
for accurate emitter position estimation. 

For ground forces, passive location finding becomes important for locating 

enemy combatants or terrorists in mountainous areas, such as those found in the 

southeast part of Turkey. Turkish land forces have been conducting search and 

destroy operations for over 30 years in this region. As a part of these operations, 

Turkish Land Forces typically first attempt to locate the terrorists who hide within 

the cave infrastructure common to the area and then destroy them. Terrorists in 

the region are known to move as small groups and to use frequency modulation 

(FM) based communication devices. If friendly forces can use a passive 
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technique to locate them this would make the search process easier. Depending 
upon the operation, an accuracy of 100-200 meters is sufficient. For this 
purpose, unmanned air vehicles (UAVs) can be used. Since typical insurgent FM 
communications are not long-duration or long distance, troops should have 
mobile direction finding (DF) devices with them or nearby. 

High clutter brings out another problem. When troops use DF devices in 
mountainous areas such as southeast Turkey, devices might receive signals 
more than once due to the reflection of the signal from the clutter. This makes the 
problem more complicated and more difficult, especially when troops want to 
locate the enemy emitter in three dimensions. 

1. Tactical Situation 

This thesis is concerned with a tactical situation in Turkey. 

Turkish armed forces has been planning and conducting military 
operations against the PKK (Partiya Karker Kurdistan (Kurdistan Worker's 
Party)). The PKK, whose main purpose is to weaken the ruling government in 
Turkey, is classified as a terrorist organization. Most nations of the world accept 
this classification. 

These operations are conducted mostly in the southeast region of the 
Turkey. That part of Turkey is very mountainous and includes many caves. 
Terrorists hide in the caves and come out whenever they want to attack. 
Operating within mountainous terrain gives them an advantage in hiding from 
government forces and helps them to move more discretely. If military forces try 
to follow them, especially at night, they escape easily because they know the 
terrain better than the military forces. Darkness covers them and mountains 
hinder the movement of the larger military forces more than the terrorists. 

These terrorists mostly come from Iraq, Iran and some other neighbor 
countries where they receive their training. They learn how to fight and use 
mines, explosives and other weapons and equipment. They cross the border in 
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small groups at night. They travel to their destinations on foot, sometimes taking 
from one to two weeks. 

The terrorists typically use hand-held radio devices for communication 
between groups and the main camp if they are close to it. Some high-ranking 
terrorists use cell phones or satellite phones if they can, but GSM companies do 
not cover most of the mountainous areas where they operate. The devices used 
between groups normally use FM modulation. They try to use these devices as 
little as possible to avoid detection and location. Often the duration of these 
communications number in the seconds or at best tens of seconds. 

FDOA cannot be used because of its long duration signal requirement. 
TDOA therefore, since it does not have the same requirement, might be the best 
solution. 

Turkish forces have been using some technology and intelligence support 
to find the terrorists, but they still manage to escape using the methods 
discussed above. 

These terrorists must be found and destroyed. A better method to find 
them is to use some sort of geo-location system targeted against their 
communications signals. Even though they try to use FM as little as possible, 
TDOA is a very good technique for this purpose. If land forces can use a tactic 
based on TDOA they will be more successful in finding, locating, and destroying 
the groups of terrorists. 

2. TDOA 

Time Difference of Arrival is one of the basic DF techniques used to locate 
RF emitters for CW signals. “TDOA takes advantage of the fact that a transmitted 
signal will arrive at different sensors at different times.” (Batson, McEachen, & 
Tummala, 2012) “A number of spatially separated sensors capture the emitted 
signal and the time differences of arrival (TDOAs) at the sensors are determined. 
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Using the TDOA’s, emitter location relative to the sensor can be calculated.” 
(Chan & Ho, 1994) 

B. MAJOR RESEARCH QUESTIONS 

This thesis will explain the following subjects. 

• What are EW and its components? 

• What is the role of DF in EW? 

This study will answer the following questions: 

• What are the advantages and disadvantages of using TDOA? 

• How effectively can ground troops use TDOA DF techniques in a 
high clutter area with noise? 

• How can we simulate TDOA technique? 

• What are the optimal uses of the TDOA technique? 

C. IMPORTANCE AND BENEFITS OF THE STUDY 

This study can be used as a guide for developing a geolocation system 
based on TDOA. It will also identify the benefits of using a TDOA-based 
geolocation system for ground forces in high clutter terrain. It will clarify the major 
EW components and their usage as an introduction to discussions of TDOA. 

D. ORGANIZATION OF THE THESIS 

This thesis is composed of five chapters. Chapter I introduces the area of 
research, major research questions, importance and benefits of the study and the 
organization of the thesis. 

Chapter II presents information about EW fundamentals. EW and its 
components are examined. The relationship between EW and 10 is explained. 
The importance of EW is explained. 
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Chapter III explains DF techniques such as directional antenna, phase 
interferometry, TDOA, and FDOA. TDOA, as a part of ES, is explained briefly. 
The details of the TDOA technique and its usage in the field are examined. 
Capabilities of TDOA are described. The developed TDOA-based location 
system for ground forces in a high clutter area are explained from Ezzat’s paper. 
The least square estimation is explained and the required recursive least square 
estimation approach is examined in detail. Sources of error are explained for DF 
and TDOA. 

Chapter IV explains the developed simulation for the system. Assumptions 
and restrictions for the simulation are listed. After running the simulation, five 
types of analyses are done to see their effects on the accuracy of the proposed 
system and the simulation. 

Chapter V is the conclusion. The results obtained from the system 
simulation are explained. Recommendations for future works are given. 
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II. ELECTRONIC WARFARE 


In this chapter, information operations and its relationship to electronic 
warfare are defined. Then electronic warfare and its components explained. 
Since the topic of the thesis is direction finding, and because direction finding is a 
subdivision of electronic warfare support (ES), ES is explained in more detail. 
These definitions are made according to Frater and Ryan’s book on EW, with 
secondary reference using U.S. EW doctrine as discussed in Joint Publication 
3.13.1. 

Another classification of EW components is based on their activities, 
which can be defined as either active or passive. Active/passive classification is 
discussed. 

A. INFORMATION OPERATIONS 

Information has a vital importance for any nation and for its existence. 
Every nation around the world tries to obtain information, use it as much as 
possible for its own good, and prevent the enemy from taking advantage of it. 
Information operations (10) is the effort done for this purpose, by human or by 
machine. 10 incorporates the actions taken to preserve the integrity of one’s own 
information system from exploitation, corruption, or disruption, while at the same 
time exploiting, corrupting and destroying the adversary’s information system 
(Adamy, 2004). 

Using U.S. doctrine as a reference, 10 coordinates and synchronizes five 
core capabilities to help the commander to reach his/her purpose in the 
battlefield. These capabilities are psychological operations, military deception, 
operations security, electronic warfare and computer network operations (JP 
3.13, 2006). In addition to the core capabilities doctrine defines supporting 
capabilities, which are information assurance (IA), physical security, physical 
attack, counterintelligence, and combat camera, (JP 3-13, 2006) and related 
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capabilities of civil military operations, public affairs and defense support to public 
diplomacy (JP 3-13, 2006). 

The world has gone to cellular and wireless. Since the Information Age 
produced a revolution in military operations, the electromagnetic spectrum is 
essential to the transmission and reception of information in the modern world, 
and therefore EW is a critical element of information operations. As 
communication and information systems become increasingly vital for military 
and civilian society, they can become targets in war for an enemy; therefore, 
they can play a significant role for offensive and defensive operations. The 
military has adopted communication and information systems and they have 
become essential for military operations. Commanders need information to be 
able to cope with the complexity of modern warfare. This reliance on information 
in turn makes them vulnerable to attack. There is an emphasis for commanders 
to attack adversary information and communication systems. Modern battlefields 
rely heavily on the use of the electromagnetic spectrum (EMS), whether for 
surveillance and target acquisition, passage of information, processing of 
information, or destruction of enemy forces (Frater, & Ryan, 2001). Figure 1 
shows all the necessary capabilities for 10 and their roles. 
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Figure 1. Information Operations Integration (From: JP 3.13) 

B. ELECTRONIC WARFARE AND ITS COMPONENTS 

Electronic Warfare is one of the five core capabilities of 10. EW dominates 
the EMS for the benefit of friendly forces. The EMS is shown in Figure 2 and 
Figure 3. EW can be defined as "... any military action involving the use of 
electromagnetic (EM) and directed energy to control the EM spectrum or to 
attack the adversary.” (JP 3.13, 2006) The means to conduct EW targeting are 
typically technological in nature, and the immediate targets are normally 
technological, but the ultimate target is the commander’s decision-making 
capability. If EW is successful, the friendly commander can make good decisions 
based on the information coming from technological devices; on the other hand, 
the adversary commander cannot make good decisions because he/she cannot 
access good and useable information as much as he/she needs to. 
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THE ELECTROMAGNETIC SPECTRUM 



Figure 2. Electromagnetic Spectrum (From JP 3.13.1) 



Figure 3. Electromagnetic Spectrum (From: National Aeronautics and Space 

Administration, Science Mission Directorate) 
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All components of EW can be applied to all kind of operations. During 
peacetime, military forces try to use EW to detect potential adversary EMS usage 
and gather intelligence; during wartime, they use EW to protect their own EMS 
usage ability and prevent adversary usage. 

EW can be divided into two parts: communications EW and non¬ 
communications EW. Communications EW is mostly concerned with 
communication sources that transmit in frequency bands between HF (3-30 
MHz), VHF (30-300 MHz), UHF (300-3000 MHz), and SHF (3000 MHz to 30 
GHz). Non-communications EW is mainly concerned with radar systems, of 
which some operate in the lower communications bands, but most are located in 
the higher microwave and millimeter wave frequencies. Additionally, research in 
directed energy weapons and technology is increasing in importance. 

As shown in Figure 4, EW has three doctrinal subdivisions. 



Figure 4. Subdivisions of EW 


1. Electronic Attack 

EA involves actions taken against personnel, equipment or facility to 
degrade their use of the EMS and combat abilities. 

Subdivisions of EA are shown in Figure 5. 
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Jamming 



Figure 5. Subdivisions of EA (From Frater & Ryan, 2001) 

2. Electronic Protection 

EP involves actions taken to protect personnel, equipment and facilities 
from any effect of friendly or adversary EW activities that degrade, neutralize or 
destroy friendly combat capabilities. 

Subdivisions of EP are shown in Figure 6. 


14 







Siting 



Figure 6. Subdivisions of EP (From Frater & Ryan, 2001) 

3. Electronic Warfare Support 

Since this thesis focuses on locating adversary emissions in the difficult 
terrain of eastern Turkey, our discussion of ES will be more extensive and 
relevant to the thesis than the other doctrinal subdivisions of EW. ES involves 
actions taken to identify, intercept and locate intentional or unintentional radiated 
electromagnetic energy. The purpose of the ES is target recognition. The main 
functions of ES are to produce intelligence, to produce steerage for EA and to 
cue surveillance and target acquisition resources. (Frater & Ryan, 2001) 

Subdivisions of ES are shown in Figure 7. 
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Search 



Figure 7. Subdivisions of ES (From Frater & Ryan, 2001) 

ES differs from, but is similar to Signals Intelligence (SIGINT). ES is used 
for immediate battlefield information and SIGINT is used for intelligence. ES 
supports near term operational applications and SIGINT supports long term 
applications. The combat information gathered by ES can be provided to 
intelligence resources in addition to being used operationally. Combat 
information does not normally require the type of deep analysis that SIGINT 
typically does. 

Previously shown in Figure 7 were the subdivisions of ES, which we will 
now discuss in more detail. 

Search: it is necessary to search for and identify the EM signal that the 
adversary uses in the EMS before it can be examined. 

The search systems act in space, time, and frequency. They have to be 
close enough to the adversary’s transmitting system to be able to detect the 
signal. They have to be actively searching or listening at the same time that the 
adversary system is transmitting. Finally, they must be listening to the same 
frequencies used by the transmitter. This frequency requirement is generally met 
with two compatible technology approaches, narrowband receivers and 
broadband receivers. Narrowband receivers can receive a single signal at a time 
and can scan a desired bandwidth sequentially in frequency; on the other hand, 
broadband receivers can monitor multiple channels at the same time. 
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Intercept: Once the signals of interest are identified through the search 
process they have to be analyzed during intercept based on their modulation, 
bandwidth, amplitude, frequency and other parameters. This process is also 
called monitoring. 

Direction Finding (DF): The location of the transmitter is determined by 
information gathered during the search and intercept process. These locations 
are likely to be an approximation rather than an exact location. DF was 
historically based on triangulation where there had to be at least three receivers 
around the emitter. 

DF systems historically employed special antennas, which defined the 
bearing towards the emitter. When the lines of bearing from each receiver are 
drawn on a map manually or automatically, the interception of the lines form a 
triangle. The smaller the triangle, the better the accuracy of the system. The 
emitter lies inside the triangle. 

Analysis: Once the signals are examined they are analyzed to define the 
adversary’s electronic warfare capabilities. The main purpose is to clarify the 
battlefield for the commander from the EMS perspective. 

ES targets an adversary’s EA, communication systems and electronic 
systems. A typical target is an adversary’s communication systems, where the 
information gathered is used for operationally actionable intelligence and 
targeting purposes. 

Because JP 3.13.1 is the shaping document of EW for the U.S. and many 
of its allies, it is a good idea to supplement Frater and Ryan’s book with 
discussions from the joint publication. According to JP 3.13.1, EA has five 
subdivisions: Electromagnetic Jamming, Electromagnetic Deception, Directed 
Energy, Anti-radiation Missiles and Expendables (e.g., chaff, flares and active 
decoys). ES has three subdivisions: Threat Warning, Collection Supporting EW 
and Direction Finding. EP has three subdivisions: Spectrum Management, EM 
Hardening and Emission Control. Subdivisions of EW do not act alone, they 
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interact with each other. The interaction among the subdivisions of EW and their 
subdivisions are shown in Figure 8. 



Figure 8. Overview of EW (From JP 3.13.1,2007) 

Electronic warfare can also be categorized by whether it is considered 
active or passive. ES tends to be passive, EA tends to be active and EP tends to 
be both passive and active. Figure 9 shows this relationship. Active activities 
require emission of detectable signals by the party conducting EW that are 
transmitted (such as in the example of the jamming of a radar). Passive activities 
do not emit signals, but rather depend upon detection of signals emitted by a 
targeted emitter. Active activities can be normally be implemented during 
peacetime only under strict limitations; on the other hand, passive activities can 


18 








be implemented during peacetime with few if any limitations (Frater & Ryan, 
2001 ). 



Figure 9. Categorization of EW based on Active and Passive (Frater & Ryan, 2001) 

In this Chapter we have discussed the importance of Electronic Warfare. 
In the next Chapter we develop the fundamentals of geolocation that will be 
important to our study. 
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III. EMITTER GEOLOCATION 


A. INTRODUCTION TO EMITTER GEOLOCATION 

Determining the location of a target emitter is one of the fundamental 
operations of EW, and can serve many useful purposes. The position of an 
emitter may indicate the position of the enemy forces. In addition, precise 
location of the target emitter enables the use of Global Positioning Systems 
(GPS) based weapons. 

For civil purposes, knowing the positions of nodes in a wireless network 
for commercial uses enables a variety of functionalities, such as emergency 
services, identification and tracking, location dependent computing, health 
monitoring and geographic routing (Xu, Ma, & Law, 2006). 

Passive location finding, in addition to active location finding techniques, 
can locate stationary and moving targets by measuring the electromagnetic 
radiation emitted by a target with the added benefit of not having to radiate 
electromagnetic energy to locate it. Passive and active location finding 
technologies play important roles in navigation, aviation, aerospace and 
electronic warfare (Yan-Ping, Feng-Xun, & Yan-Qui, 2010). 

The purpose of direction finding (DF) is to estimate or fix the position of 
selected emitters. This position is not certain because all the measurements 
include some sort of error, and the entire system includes the noise that is found 
in all communications systems. 

Several techniques can be used to calculate the position of a target 
emitter. These techniques are based on different types of information acquired 
from the received signal to calculate a position fix (PF). 

The azimuth angle of arrival of a signal, or so called line of bearing (LOB), 

is the most commonly used technique for calculating a PF. Two or more LOBs 

are used to determine a position in two dimensions (2D). These LOBs are 

assumed to be measured at the same time on the same target. These LOBs may 
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intersect as illustrated in Figure 10. This technique is called triangulation (Poisel, 
2005). 

Since two LOBs intersect at a point, the information to fix the position of 
the emitter is not accurate due to measurement and propagation errors. For 
triangulation, at least three receivers are located on a baseline as illustrated in 
Figure 10. Each DF receiver has special antennas. These antennas are used to 
measure the bearings. These bearings are plotted on a map either manually or 
automatically. Intersection of these bearings (lines) forms a triangle and the 
possible location of the target is calculated to be at the middle of the triangle. The 
size of the triangle depends on the accuracy of measurements. The smaller the 
triangle, the better the accuracy (Frater & Ryan, 2001). 



Receiver 2 


Figure 10. Intersection of measured LOBs 

Another emitter location technique is to measure the time of arrival (TOA) 

of a signal at several sensors. The TOAs can be used to calculate the position 

directly, but typically they are sent to a central node where the time difference of 
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arrival (TDOA) is calculated from every pair of TOA. Then the range differences 
between sensors and the target are calculated. These range differences are 
related to the TDOAs by the speed of propagation in the medium. In the air, this 
is assumed to be the speed of light. 

The TDOA technique generates quadratic lines of position (LOP). All the 
LOPs are subject to measurement errors and noise. The intersection of these 
LOPs is used to define the emitter position (Poisel, 2005). 

The next step of DF is geolocation. Geolocation is closely related to DF 
and triangulation, but it is more realistic and distinguished from DF by 
determining a meaningful location rather than just a set of geographic 
coordinates. 

“Emitter geolocation has two components. One is measurement, or choice 
of sensors, and the other is estimation/information fusion, or processing of 
measurements provided by the sensors.” (Musicki & Koch, 2008) 

Geolocation of a source has a wide variety of applications, such as 
location of radar sites. Localization of a interference source in satellite 
communication systems is another example of geolocation (Sathyan, Kriburajan, 
& Sinha, 2004). 

Geolocation is based on techniques that rely on frequency, time, or spatial 
information, or a combination of these. Common methodologies use angle of 
arrival (AOA), TOA, TDOA, or differential doppler (DD), also called FDOA (Ho & 
Chan, 1993) (Loomis, 2007). 

B. BEARING ESTIMATION 

Several techniques that can be used to determine the LOPs are discussed 
in this section. The phase or time difference can be measured between the 
signals. The amplitude difference between two signals can also be measured. 
Frequency differences measured can be used to determine the bearing, if one or 
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more of the receivers are moving relative to the other or to the target (Poisel, 
2005). 

LOP systems do not operate at the frequencies of the signals. The 
frequency is typically converted to an Intermediate Frequency (IF) and the 
phase/time measurements are made on this converted signal (Poisel, 2005). 

1. Circular Antenna Array 

One of the common types of antenna arrays for bearing determination is a 
circular array. An example of a four-element array is illustrated in Figure 11. 
Other forms of this antenna may include more or fewer elements. A sense 
antenna is included with a circular array. The sense antenna is used to remove 
ambiguities (Joong, Chul-Gu, & Gyu, 2004). 

Let R be the radius of the circular array. R is the length of an “arm” of the 
array measured from the center to any one of the antenna elements. The 2R/A is 
referred to as the aperture of a circular array (Poisel, 2005). 



Figure 11. Four Element Circular Antenna Array (From Poisel, 2005) 
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An incoming signal s(t) and its accompanied E field with a magnitude E is 
illustrated in Figure 12. The vertical component of this E field is the only portion 
that affects a vertically oriented antenna. Thus, E vert is given by Ecosa and the 
corresponding signal amplitude must be adjusted by this factor (Poisel, 2005). 



2. Interferometry 

One of the techniques for measuring the AOA of a coming signal, and 
determining its LOP, is interferometry. In interferometry the phase difference or 
time difference between two antennas is measured directly. The distinction 
between this technique and other techniques lies in what is measured. 
Interferometry requires highly accurate measurements. In order to achieve this 
performance accuracy it is necessary to space at least two of the antennas such 
that the range of possible phase difference can exceed 2tt radsian (Vaccaro, 
1993). 
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Phase comparison DF systems consist of several antenna elements which 
are arranged in a particular geometric configuration. The number of the elements 
and the arrangement depends on the DOAs of interest and the method used to 
process the signals. The minimum number of elements is two. The determination 
of DOA is performed by direct phase comparison of the received signals from the 
different antenna elements. 

There are two types of interferometers. The first type measures the phase 
difference between the two antennas and calculates the AOA from that 
measurement. These interferometers are called phase interferometers. The 
second type measures the time of arrival differences to the two antennas and 
calculates the AOA from these measurements. These types of interferometers 
are called active time interferometers. For this second type of interferometer 
there must be some sort of time mark to be able to measure the time difference. 
If there is no time mark then the signal reaching these two antennas must be 
correlated. Radar pulses provide a convenient time mark in their leading edge. 

Interferometer systems operate well over limited frequency ranges. They 
provide the capability of receiving signals with good accuracy (Lipsky, 1987). 

C. QUADRATIC POSITION FIXING TECHNIQUES 

Techniques for determining a position fix based on Time of Arrival (TOA), 
Time Difference of Arrival (TDOA) and Frequency Difference of Arrival (FDOA) 
(also known as Differential Doppler or Differential Frequency) techniques are 
explained in this section. 

When TOA or TDOA is measured at two or more widely separated 
receivers, quadratic LOPs are determined. The intersection of these LOP curves 
is taken as the estimated location of the emitter. 

The receivers and the emitter might be stationary for TDOA and TOA but 
for FDOA either the receiver or the emitter must be moving in order to produce a 
frequency difference induced by movement, necessary for measurement. 
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The main advantages of using these techniques are: 

• Most of the time only one antenna is used per receiver instead of 
two or more antennas as in interferometric techniques. 

• Higher precision and more accuracy can be obtained with these 
techniques (Poisel. 2005). 

On the other hand, there are two main disadvantages: 

• Preprocessed data samples are required for non-pulsed 
modulation, like Continuous Wave (CW) signals such as using 
Frequency Modulation or Amplitude Modulation (Poisel, 2005). For 
pulsed conventional radar signals it is easy to define the beginning 
of the signal, but for CW signals more complicated methods like 
correlation should be used. 

• Also, it is difficult to measure frequency of pulse-type signals to the 
level of accuracy needed to do FDOA because the frequency 
resolution is equal to 1/7", where Tis the pulse duration. 

This following section begins with a presentation of the TDOA technique, 
followed by a discussion of FDOA, and concludes with a discussion of TOA. 

1. Time Difference of Arrival (TDOA) 

The TDOA involves the measurement of the TOA of the received signal. 
The TDOA technique needs two or more geographically separated sensors 
synchronized with each other to be able to find the location of the emitter. If only 
two receivers are used, this will normally result in an ambiguous solution of two 
locations, and a third receiver is necessary to resolve this ambiguity. These 
receivers can be either on the ground or on an airborne platform. 

The geometry of TDOA is shown in Figure 13. Only two receivers are 
shown for simplicity and will be used for calculations. This geometry is for 2D but 
if the receivers or the emitter is elevated it becomes 3D. For both the scenarios, r 
is the slant range where o is the range between receiver 1 and emitter and r 2 is 
the range between receiver 2 and emitter. 
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Figure 13. Geometry of TDOA 


The geometry involves an emitter at position x = (x T ,y T ) and two receivers 

at positions (-^,0) and (x 2 ,0), At time t = 0, a measurement of the TDOA is 

made between the arrival of the same pulse from the emitter arriving at the two 
receivers. 

The distances between receivers and the emitter can be calculated as 
follows: 


r, = ct t i = 1,2 


(3-1) 


where; 

c is the speed of light, 

tj is the time between when the signal leaves the emitter and 
when it arrives the receiver. 
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TDOA is the time difference between when the signal arrives at one 
receiving site and at the other, and represented as r 

r = 4-*i=-= (3.2) 

c c c 

r,.J(x T - Xf + y T 2 '=1.2 (3.3) 

re = 5= ^(x r -x,) 2 + y r 2 -^(x r -x 2 ) 2 + y r 2 (3.4) 

Squaring both sides of the previous equation yields: (the receivers are at 
the same distance from the origin) 

5 2 = 2x r 2 + 2y/ + 2s 2 - 2,|(x r + s) 2 + y r 2 ] [(x T - s) 2 + y r 2 ] (3.5) 

where 

X, = -s & x 2 = s 

With some more algebra, this expression becomes 

- 4s 2 S 2 = (4S 2 - 16 s 2 )x t 2 + 4 8 2 y T 2 (3.6) 

After further manipulation, this becomes: 
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(3.7) 


x T 


Yt 


8 2 !4 (4s -8 )/4 


This is familiar equation for a hyperbola, which has x-axis intercept of (0, 
0) and is asymptotic to the lines (Loomis, 2007) 


y = ± 


^4s 2 - 8 2 
8 


x 


(3-8) 


The curve defined by this expression is a hyperbola. Several of the 
hyperbolic curves are shown in Figure 14. 



Figure 14. Hyperbolic TDOA Curves (From Loomis, 2007) 
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It is clear that two receivers cannot find the location of the emitter since 
these hyperbolas covers a wide variety of locations. At least three receivers are 
necessary to find the geolocation of the emitter in two dimensions or on the 
earth’s surface. 

The following derivation of the solution for the location of the emitter as a 
mathematical model of TDOA is quoted from Poisel’s book, Introduction to 
Communication Electronic Warfare Systems, Chapter 12. 

Suppose there are S receivers available to receive and compute the 
location of the emitter, then the Equation (3.1) becomes for all pairs of sensors 
(Poisel, 2002). 


d i -d j = c(t i -t j ) = ct j 


/,y = 0,1,...S-1 (3.9) 


According to Poisel’s approach, let’s assume that all of the arrival times 
are compared with the arrival time at a sensor located at coordinates (0,0) as 
shown in Figure 15. 
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Figure 15. Sensor Grid and target in Two Dimensions 


The time differences of arbitrary (i,j) are not used, just (/,0) is used for all 
/.The Equation (3.9) in this case reduces to 


d, = IK r T 1 1 = c(t,, 0 + ^ 


(3.10) 


When putting the locations of the emitter and the receivers in to Equation 
(3.10), it becomes 


[*/ y, z i] 


X-r 


y T 


+ 


ct ifl yjx T 2 + y T 2 + Z T 2 


(3.11) 


=4 c v4 (x ' 2+y ' 2+z ' 2) 
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Putting Equation (3.11) into matrix form 


PiX T +cy x T i=“Cf /i0 


2 1 ll 
+ - P; 
2 11 


(3.12) 


where 


p, is the position vector of receiver / 
Expanding this result for all / = 1 ...S receivers yields 


Px T +c||x T || = d 


(3.13) 


where 


A *i /. z i A 


* s i y s i z 


s-1 J 


X T = 


x T 




1,0 


t = 


s-1,0 


c = ct 
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d = ^ diag{PP T - c 2 tt r ) 


Let a = ||x T || and Q = (PP r ) 1 , then the expression becomes 


(c r Qc -1) a 2 - 2d r Qca + d r Qd = 0 


(3.14) 


which is a quadratic equation that can easily be solved for a which is the 
range of the target from the origin. Substituting this range back into Equation 
(3.13) will solve the problem for the target location. 


Px T + ac = d 

x T =P ^(d-ac) 


(3.15) 


If S>4, this becomes an over-determined system of equations. Instead of 
the inverse, then, the pseudo inverse can be used and is shown in Equation 
(3.16) 


pt =(P 7 P)- 1 p r 


(3.16) 


The pseudo inverse solves for the emitter position in the minimum least 
squares error sense (Poisel, 2002). It is very important for this thesis and for 
many of TDOA solutions. The pseudo inverse is used in Chapter 4 for simulation 
which is a developed version of Ezzat’s Closed-Form Geolocation Solution. 
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2. Frequency Difference of Arrival-Differential Doppler (FDOA-DD) 

The signal emitted by a target of interest which is moving produces an 
effect called Doppler shift. The Doppler shift is related to the direction of the 
movement relative to the receiver and shows itself as a frequency difference. The 
movement can be at the target or at the receiver. Each of these movements 
produces the same frequency difference effect. When the frequency difference 
between the target and the receivers of two (the more the receivers the better the 
calculation) is measured, these measurements can be used to calculate the 
geolocation of the target. This type of geolocation technique is called frequency 
difference of arrival (FDOA) or differential Doppler (DD). 

Since this thesis is interested in stationary receivers and FDOA technique 
requires a moving receiver or a moving target, FDOA is not examined in detail. 

D. CLOSED-FORM SOLUTION OF HYPERBOLIC GEOLOCATION 

EQUATIONS 

3D geolocation is important for the Turkish Army because the terrain 
where most of the anti-terrorist operations take place is very mountainous, 
especially the southeast part of Turkey. It is not possible to locate or position the 
receivers at the same elevation with the transmitter. It is inevitable to be 
presented with a 3D problem. Ezzat’s approach gives a unique solution to a 3D 
problem from the TDOA perspective. Because of that reason, Ezzat’s approach 
is used for developing and analyzing the scenarios in Chapter 4. 

The following approach is paraphrased from Ezzat’s article. 

The approach that is explained here does not depend upon range data. 
Range data is derived by multiplying the time that the signal propagated in air 
and the speed of the signal, the speed of light. It is necessary to know the time of 
transmission and the time of arrival to get the amount of time travelled from the 
transmitter to the receiver. For the most part, it is impossible to know the time of 
transmission, which can be defined as f 0 . Ezzat indicate in his paper that this is 

the only technique that can be used without the range data. He mentions that 
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without t 0 , it is possible to find the location of the emitter. He also indicates that 
the previous solutions are good for a noise free environment. 

The closed-form solution presented does not require the calculation of 
range data and does not depend on the availability of any information other than 
the times of arrival. 

The basic form of time of arrival equation is as follows 


''=^t 


(3.17) 


where 

f, is the time of arrival at receiver /, 

D, is the distance between the emitter and the receiver, 
c is the speed of light. 
t 0 is the time of the transmission 

Two or more receivers are needed to be able to calculate the TDOA as 
mentioned before. When this condition is satisfied, it is possible to eliminate t 0 
from any pair of two equations, which results in 


= (3. 18) 

This is the TDOA equation. This equation yields a 3D hyperboloid as 
mentioned before. When the emitter coordinates, x 0 , y 0 , and z 0 , are plugged 

into this equation,(3.18) can be written as 
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>/(*2 ~ x o) 2 +(y 2 -y 0 ) 2 +( z 2 - z o ) 2 - >/(*i - * 0 ) 2 +(/i - y 0 ) 2 +( z i - z o ) 2 


= c(f 2 -0 


(3.19) 


where 


^, y 1; z, and x 2 , y 2 , z 2 are the coordinates of the receiving 
antennas 1 and 2. 

Having three more receiving antennas yields three additional times of 
arrival, which produce an additional three equations like Equation (3.19), it is 
possible to solve for the emitter coordinates x 0 , y 0 , z 0 

There is an effect of path delay on (3.18) and hence(3.19). Ezzat says that 
a propagation mode between any two points in which the path of the signal is not 
a straight line will be mathematically equivalent to propagation along a straight 
line but with a velocity that is less than c, as the time of arrival is important. The 
TDOA equation for the case in which the path of the signal is nonlinear, it will be 
written as follows after adding the path delay: 




D, 


a 2 c a^c 


1 ( 4 _9i) 

C a 2 


(3.20) 


where 

a 1 and a 2 are path delay coefficients which are less than or equal 

to 1. 

a is one when there is no path effect on the propagating signal like in the 
case where the signal propagates through air. For this thesis, all the path delay 
coefficients are assumed to be one. 
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The solution that is used for this thesis is derived from Ezzat’s article. 
First, I will explain his approach. Then, because he did not discuss the effects of 
noise, I will add noise to his approach. 

These are the steps to transform the hyperbolic equations into a set of 
vector equations. In Equation (3.18), we first note that a distance d j can be 
written as the norm of a vector. 


cHPi-Po 


(3.21) 


where 

Pi =(x l ,y i ,z i ) 

and 


is the position vector of the receiving antenna 


Po = (Wo’ z o) is the position vector of the emitter. 


Equation (3.20) represents the difference between the two. We modify 
Equation (3.20) and write it as a difference of squares 


1 rJ ^ ri ^ 
(t 2 -t 0 ) 2 -(t,-t 0 ) 2 = M\-%) 


a n 


a , 


(3.22) 


From Equation (3.21) and Equation (3.22) we have 
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(3.23) 


= c 2 (4 2 + fo 2 -2M, - f. 2 - 'o 2 + 2Vo) 

a 2 v 7 

= c 2 (t 2 2 -t?)-2t 0 c 2 {t 2 -ti) 

Pi - p 0 1 2 can be written as 


Pi-Pof = |Pi| 2 + |Po| 2 - 2Pi T Po 


(3.24) 


where 


represents the transpose of vector Pj. 

If we put Equation (3.24) into Equation (3.23) it gives 


p 2 r+ip 0 r ■ 2 p 2 t po iPi r+ip 0 r - 2 pi t p 0 


a 0 


a. 


c 2 (t 2 2 -t 2 )-2t 0 c 2 (t 2 -t,) (3.25) 


The two coefficients a 1 and a 2 are both very close to one, and as a result 

the difference |p 2 1 2 /a 2 2 - |p 1 1 2 la 2 is negligible by comparison with the other 
quantities. Then Equation (3.25) becomes 


P 2 Pi 


a n 


a, 


2 Po ~ ^V) = c " (4 2 - 0) ~ 2f 0 <c 2 (4 - 0 

a 2 a, v ’ 


(3.26) 
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In Equation (3.26) there are two unknowns, t 0 and p 0 . To solve the 
coordinates of the receiver (p 0 ) we need three linearly independent equations 

like Equation (3.51). Equation (3.26) shows the relationship between receivers 1 
and 2. The other two linearly independent equations can be between 1 and 3 and 
1 and 4. This pairing of receivers makes the equations linearly independent. As 
can be seen for three equations at least four receivers are required. For the 
closed-form approach, it is necessary to get rid of t 0 . To do this these linearly 

independent equations are going to be coupled to subtract from each other. For 
that purpose this solution needs five receivers to work. 

Equation (3.26) is divided by {t 2 -t,) , then it becomes 


1 

( \ I 2 

lp 2 l 

1 I 2 ^ 

|Pi| 
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CM 
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(n T 

P 2 
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('2-4) 

a 2 2 
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a. 2 

1 J 

1 (4-0 1 

V a 2 

«i 2 v 


= c 2 {t 2 + t,)-2t 0 c 2 


(3.27) 


By similar steps Equation (3.27) can be expressed for receiver pairs 3-1 
and 4-1 
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( 1 I 2 
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1 

( _ T 

P 3 
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Equation (3.28) is for 3-1 receiver pair. 

1 

(\ I 2 1 |2 X 

|P 4 | |Pl| 

2p 0 

f _ T n T\ 

P 4 Pi 


a 4 a , 2 

v 4 1 J 

('4-0 

2 2 
a i ) 


c 2 {t 3 + t,)-2t 0 c 2 


c 2 {t 4 + t,)-2t 0 c 2 


(3.28) 


(3.29) 
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Equation (3.29) is for 4-1 receiver pair. 

The next step is to eliminate t 0 from Equations (3.27), (3.28) and (3.29). 
(3.27) and (3.28) are merged in the components of p 0 to get rid of t 0 . This results 
in 
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f | |2 

|P 2 | 

1 I 2 Y 

|Pl| 
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( \ |2 
|P 3 | 

i i 2 A 
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(3.30) 
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Similarly Equations (3.27) and (3.29) yield 
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(3.31) 


The last two equations are linearly independent in the components of p 0 . 

As mentioned before at least three independent equations are required. For this 
purpose there must be another receiver: receiver 5. If we follow the same steps 
that we followed to reach Equations (3.30) and (3.31), we can have the following 
equation for receiver-emitter pair 5-1. 
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(3.32) 


Equations (3.30), (3.31) and (3.32) can be written in alternative algebraic 

form 


3,1*0 ^12/0 ^13^0 

^21*0 ^ 22 X 0 ^23 Z 0 — ^2 
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(3.33) 
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a. 


(3.35) 


The other constants can be derived from the same Equations (3.31) and 

(3.32). 

Equation (3.33) can be expressed as 


AX = B 


(3.36) 


and can be solved as in least square estimation way as shown in Section 
E. 


X = (A T * A)' 1 A T * B 


(3.37) 


E. LEAST-SQUARE ESTIMATION 

For the least-square estimation, we use Barkat’s book, Signal Detection 
and Estimation, 2005. 

In the least-square estimation, the criterion is to minimize the squared 
difference between the given data (signal plus noise) and the assumed signal 
data. 

This development applies to a linearized version of non-linear equations 
relating an /^-dimensional measurement vector and the /W-dimensional vector to 
be estimated. This linearized matrix equation represents the difference between 
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the actual measurement vector and the measurement that would be obtained if 
the vector to be estimated has the estimated value. 

Suppose we want to estimate M parameters, denoting the /W-dimensional 
vector 0, from the K measurements, denoting the K'-dimensional vector y with K 
> M. The relation between the parameters 0 and the observed data y is given 
by the linear model 

Y = H0 + N 

where 

h is a known (K x M) matrix 

n is the unknown (K x 1) error vector that occurs in the 
measurement of 0. 

The least-square estimator of 0 chooses the values that make X = H0 
closest to the observed data Y- Hence, we minimize 


J(0) = £(Y k - X k f = ( Y - H0) r (Y -H0) = YY T - 2Y T H0 + 0 T H T H0 (3.38) 


k =1 


Note that Y T -H0 is a scalar. Taking the first-order partial derivative of the 
cost function J(9) with respect to 0 and setting it equal to zero, we obtain the set 
of linear equations 


= -2H t Y + 2H T H0 = 0 
dQ 


(3.39) 


and LSE can be found to be 


44 



0 |S =(H T H)' 1 H T Y 


(3.40) 


We observe that the error in the estimator 0 |S is a linear function of the 
measured errors n, since 


Me-0, s = e- (h t h y h t y = 0 -(h t h) _1 h t [ho+n] = -<h t h) _1 h t n ( 3 . 41 ) 

1. Recursive Least-Square Estimator 

In real time estimation problems, it is necessary to write the estimator 0 in 

a recursive form for better efficiency. Consider a situation where an estimate 0 is 
determined based on some data Y K . If new data Y K+1 are to be processed after 

having determined an estimate based on the data Y K , it is best to use the old 
solution along with the new data to determine the new least-square estimator 
(Barkat, 2005). “It is clear that discarding the estimate based on the data Y K and 

restarting the computation for a solution is inefficient. This procedure of 
determining the least-square estimate from an estimate based on YK and the 
new data Y K+1 is referred to as sequential least-square estimation, or more 
commonly recursive least-square (RLS) estimation.” (Barkat, 2005) 

Consider the problem of estimating 9 from the data vectors z M given by 
the linear model 


z M =ly? + u M 


(3.42) 


where 
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(3.43) 


z =[Y Y Y V 

is an (MK+1) collections of vectors Y,,Y 2 ,...,Y M since each vector 
Y k , k = 1, 2, Mis a (K+1) vector 

u „ = [N,N 2 ...N u Y (3.44) 

is an (MK+1) error vector, and 

h M =[^i h 2 ...h M f (3.45) 

is an (MK*n) mapping matrix relating Z M to the (n*1) parameter vector 0to 
be estimated. 

It can be shown that the RLS estimator is given by 

- ®m-i + V m [u m (3.46) 

where 

V M = C yu h M T R MM 1 (3.47) 

c is the error covariance matrix given by 
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(3.48) 


'uu 


R, R 




12 


R 


1M 


A, S 


R 


22 


R 


2/W 


R 


1/W 


R. 


2/W 


... R, 


MM 


In three dimensions, when the covariances are fixed this matrix is 


Cuu - 


Pxy a x a y 

Pxz°x°z 


Pxy®x^y 
PyzVyVz 


Pxz°x a z 

PyzVyVz 


where 

a x is the variance of x 
<r y is the variance of y 
<j z is the variance of z 

p xy is the correlation coefficient between x and y 

p xz is the correlation coefficient between x and z 

p yz is the correlation coefficient between y and z 

If it is assumed that x, y, and z are uncorrelated and 
p xy ,p xz ,p yz = 0 . In this case Equation (3.49) becomes (Poisel, 2005) 


C 


uu 


er x z 0 0 

0 < 0 

0 0 a z 2 


(3.49) 


so that 


(3.50) 
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F. 


SOURCES OF ERRORS AND MEASUREMENT ACCURACY 


1. Sources of Errors 

Direction Finding is subject to number of errors. These error sources are 
listed and described below. 

a. Equipment Error 

Modern DF equipment gives bearing with accuracy of ±2°. Hand¬ 
held tactical units have less accuracy of ±10° (Frater & Ryan, 2001). 

b. Short-baseline Error 

If the angle between bearing lines is less than 45°, the triangle of 
error for triangulation or confidence ellipse for hyperbolic geolocation systems 
like TDOA becomes significantly larger (Frater & Ryan, 2001). 

c. Co-channel Interference 

Most tactical DF systems cannot identify the difference between 
multiple received signals. When there is significant co-channel interference these 
DF systems tend to give erroneous bearing information (Frater & Ryan, 2001). 

d. Adjacent channel Interference 

“Strong signals in a channel adjacent to the one being DF’ed can 
lead to an erroneous bearing” (Frater & Ryan, 2001). 

e. Multipath Error 

In multipath error, two or more signals arrive at the receiver. These 
received signals originate from the same source but travel in a different path due 
to natural or man-made obstacles and reach the receiver at different times due to 
the difference in the distance traveled over the separate paths. 
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f. Night Effect 

Night effect is a special case of multipath effect that occurs when 
sky-wave propagation occurs at night but not during the day. This type of error 
occurs at long distance, typically over HF communication channels. 

g. Coastal Refraction 

Surface wave propagation that crosses a coastline at an angle 
other than a right angle is subject to bending caused by refraction. This may lead 
to wrong bearing determination. Coastal refraction is usually significant at 
frequencies below 10MHz (Frater & Ryan, 2001). 

h. Thunderstorms 

Thunderstorms can lead to a wrong bearing that points towards the 
thunderstorm rather than the originating transmitter. 

i. Rain 

Heavy rain may reduce received signal levels in SHF and higher 
bands. This reduces the range of DF systems (Frater & Ryan, 2001). 

2. Cross-Correlation TDOA Estimation Technique 

The TDOA position fixing technique includes two phases. The first phase 
is the estimation of the TDOAs of the signal from a source, between pairs of 
receivers through the use of time delay estimation techniques. In the second 
phase, the estimated TDOAs are transformed into range difference 
measurements between stations, resulting in a set of hyperbolic equations 
(Aatique, 1997). 

The previous sections related to TDOA calculations are for the second 

phase. 

There are two general methods for estimating the TDOAs. The first one is 
to subtract TOA measurements from two stations to produce a relative TDOA. 
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The second one is to employ a cross-correlation technique, in which the received 
signal in one station is correlated with the received signal at another station. 

Because it is very difficult to know the timing reference on the source to be 
located, and because the signals of interest for this thesis are CW, the cross 
correlation technique is commonly used to estimate the TDOAs. The basic timing 
requirement for this technique is to synchronize the receivers. This requirement 
is relatively easy compared to the need to know the originating transmission time 
of the signal. Therefore, we will focus on the cross correlation TDOA estimation 
technique. 

Signal, s(t), emanating from a remote source through a channel with 
noise, the general model for the time-delay estimation between received signals 
at two base stations, x^and x 2 (t ), is given by (Knapp & Carter, 1976) 


x 1 (0 = s (0 + a ? 1 (0 
x 2 (t) = As(t-r) + n 2 (t) 


(3.51) 


where 

n,(t) and n 2 (t) are noises 

t is the TDOA between the receivers 

A is the amplitude ratio for scaling the signal 

This model assumes that s(t), n^(t) and n 2 (t) are real and jointly 
stationary random processes. The signal s(t ) is assumed to be un correlated 
with noise n,(t) and n 2 (t). 

The cross correlation of this two received signal is given by (Knapp et al, 

1972) 
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^(r) = E[x,(Ox 2 (f-r)] 


(3.52) 


where e represents the expectation. Equation 3.51 can also be expressed as 
(Aatique, 1997) 


= FP^t) = \^x,{t)x 2 {t-T)dt 


(3.53) 


Because the observation time cannot be infinite and can be estimated 
from a finite observation time, an estimate of the cross-correlation is given by 


^(r) = yj o r x 1 (0x 2 (f-r)df 


(3.54) 


where t represents the observation interval. The time delay causing the 
cross correlation peak is an estimate of TDOA, r . 

The cross correlation technique is affected by many errors, which can be 
considered as a Gaussian distribution. Standard deviation for that Gaussian 
distribution is explain in Section 3. 

3. Standard Deviation 

When errors producing ambiguity are taken into consideration, TDOA 
isochrones are no longer clear functions; they form regions or areas where the 
target should be, as illustrated in Figure 16. 

All of the TDOA methods are subject to errors in measurement. The noise 
and the measurement error are the two primary error sources (Poisel, 2005). For 
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the purpose of this thesis, errors discussed from this point forward will be limited 
to noise-induced errors, and they are discussed in terms of standard deviation. 



Figure 16. TDOA Isochrones with Errors 

The represented TDOA hyperbolic isochrones can be shown in detail as in 
Figure 17. The standard deviation increases close to the edges of the hyperbolic 
isochrones. 



Figure 17. TDOA Isochrones with Standard Deviation 
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The standard deviation discussed in this chapter result from cross¬ 
correlation measurements. 

The Cramer-Rao bound on parameter estimation is a frequently used 
measure on how well such a parameter can be measured. The Caramer-Rao 
bound for estimating the time of arrival of a signal at a receiver is given by (Stein, 
1981 ). 


1 1 

IWr 


( 3 . 55 ) 


where 


B is noise bandwidth of receivers 

T is integration time, which must be less than or equal to signal 

duration 

y is the effective input SNR at two sensor sites 

RMS radian frequency is given by p which is the measure of the 
bandwidth of the signal an given by 


p = 2x 


OT'iTtf 2 


( 3 . 56 ) 


where 


S(f) is the spectrum of the signal 


Variable y is a composition SNR at the two sensors and is given by 
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(3.57) 
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where 


y i and y. are the Signal to Noise Ratio (SNR) at two receivers. 
For low SNR, standard deviation is given by 
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(3.58) 


where 

T is the integration time 

8 is the SNR 

W is the bandwidth and is given by W - f 2 - 

f 0 is the center frequency 

For high SNR, standard deviation is given by 


cr. 


1 1 1 

4x 2 T 8 Jf* _ ft 


(3.59) 


The standard deviation function for low SNR is illustrated in Figure 18 and 
Figure 19 for short integration time and high integration time, and in Figure 20 
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and Figure 21 for high SNR for the 25MHz bandwidth (Poisel, 2005). The dotted 
lines represent the function in Equation (3.55) and the solid lines represent the 
function in Equation (3.59). 



Figure 18. TDOA Standard Deviation for Low SNR, W=25 MHz and Long Initegration 

Time (From Poisel, 2005) 
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Figure 19. 


Figure 20. 


TDOA Standard Deviation for Low SNR, W=25 MHz and Short Integration 

Time (From Poisel, 2005) 



SNR (dB) 

TDOA Standard Deviation for High SNR, W=25 MHz and Long Integration 

Time (From Poisel, 2005) 
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Figure 21. TDOA Standard Deviation for High SNR, W=25 MHz and Short Integration 

Time (From Poisel, 2005) 

4. TDOA Dilution of Precision 

TDOA measurement suffers from another type of error, caused by the 
long ranges from the sensor baseline. Consider the Figure 22 where the target is 
very far from the baseline between the sensors. The hyperbolic LOPs are nearly 
parallel to each other and make the measurement more difficult, thereby making 
the system vulnerable to any measurement or noise error. This is called 
geolocation dilution of precision (GDOP)( Poisel, 2005). 
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Figure 22. Geometry of GDOP Effect on TDOA (Poise!, 2005) 


The GDOP effect is illustrated in Figure 23 where the effect of the distance 
to the error can be seen easily. N represents the number of receivers in the 
system. 


58 




Figure 23. TDOA GDOP from 4 to 10 elements (From Bard & Ham, 1999) 

5. Effects of Movement on TDOA 

When there is a motion between the receivers or the target the error of 
TDOA measurements increases. This is illustrated by Chan and Ho, who called 
this effect scale difference of arrival (SDOA). Figure 24 shows the relationship 
between the velocity and the error caused by this velocity. As can be seen, when 
the velocity increases the mean square error also increases. 
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Figure 24. SDOA Caused by the Relative Motion (From Chan & Flo,2003) 

In this Chapter, we discussed the fundamentals of emitter geolocation, 
bearing estimation, quadratic position finding methods, closed-form geolocation 
of emitter based on TDOA, least-square estimation method and sources of error 
for geolocation. In the next chapter, we developed a Matlab 1 simulation based on 
Ezzat’s closed-form geolocation technique and made four type of analysis; the 
effect of the distance and noise on accuracy, comparison between stationary 
receivers and moving receivers based on their effects on accuracy, the optimum 
distribution of stationary receivers for best accuracy. 


1 Matlab is a registered trademark of The Mathworks, Inc. 
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IV. SIMULATION OF THE CLOSED-FORM GEOLOCATION 

TECHNIQUE 


The simulation described in this chapter focuses on the use of time 
difference of arrival (TDOA) employing Ezzat’s proposed closed-form solution for 
the emitter position in three dimensions (Ezzat, 2007). We will analyze its 
capabilities in terms of the effect of the number of experiments on accuracy, the 
effect of a noisy environment on accuracy, the effect of distance between emitter 
and receivers on accuracy, and the optimum positions of the receivers for best 
accuracy. Finally, we will make a comparison between moving receivers and 
stationary receivers. 

The assumptions of and the restrictions to the simulation are given in 
Section A. The development of the simulation and the followed procedures are 
explained in Section B. Simulation results are given in Section C. 

Ezzat’s closed-form geolocation technique has been used to develop this 
simulation. This technique is explained in Section.D in detail. Matlab has been 
used for simulation, and Matlab code can be found in Appendix A. 

In short, simulation is done to accomplish these analyses: 

• The effect of number of experiments on the accuracy of the 
estimated emitter position, 

• The effect of the distance between emitter and receivers on the 
accuracy, 

• The effect of the noise on the accuracy, 

• Comparison between moving receivers and the stationary 
receivers, 

• The optimum geographical distribution of the receivers for best 
accuracy. 

Different codes are developed for each analysis and can be found in the 
Appendices. 
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A. 


ASSUMPTIONS AND RESTRICTIONS 


The following assumptions and restrictions apply to the development of 
the technique. 

1. Assumptions 

1) The sensor locations are known exactly. 

2) The transmitter location is fixed during the period of DF fixing. 

3) The sensors are properly located and operated. 

4) The emitter is within the range of the receivers. 

5) TDOA is measured using some sort of technique and known. 

6) The error in the noise is distributed as a Normal distribution with 
zero mean and known standard deviation. 

7) The noise for each pair is independent. 

8) There is no multipath effect. 

9) The emitter location estimate error is distributed as a normal 
distribution. 

Assumption #1 is reasonable based on the fact that any such position 
errors can be added to the emitter estimate uncertainty, if they are significant. 

Assumption #2 is necessary to the analyses of the systems considered in 
this thesis, and it is reasonable over the period required to obtain a single fix. 
Further, any change in emitter position is not significant in comparison to the 
speed of signal transmission (speed of light) and time to calculate solutions. 

Assumption #3 is reasonable in the absence of contradictory information. 

Assumption #4 is reasonable based on the fact that if the emitter is out of 
the range, then the system will not function. 

Assumption #5 is reasonable because the measurement of TDOA is out of 
the interest area of this thesis. 
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Assumption #6 is reasonable based on the nature of the noise in every 
receiver system. 

Assumption #7 is reasonable based on the fact that every receiver pair 
has different SNR level, integration time and bandwidth. 

Assumption #8 is reasonable based on the fact that multipath effect is out 
of the interest area of this thesis. 

Assumption #9 is usual when considering measurements which are 
subject to random measurement error. There are biases in the measurements 
from navigation errors, errors in the calibration tables, interference, etc.; these 
biases may be removed. In the absence of specific knowledge about these errors 
the normal assumption is reasonable. Biases will not be treated in this thesis. 

2. Restrictions 

In addition to the assumptions discussed in this section, this thesis does 
not consider the following effects; 

1) Geographic transformation, map projection effects, and grid 
reference system conversions. 

2) Propagation effects. 

3) Susceptibility to deception. 

4) Special problems associated with low-probability-of-intercept 
emitters (low SNR, spread-spectrum, time-frequency diversity, frequency agility, 
etc.). 

5) Numerical computation and normal truncation effects. 

6) Equipment errors. 

7) Interference. 

8) Night effect. 

9) Coastal refraction. 
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10) Thunderstorms and precipitation. 


B. DEVELOPMENT OF THE SIMULATION 

Ezzat’s technique and Matlab code solve for the position of the emitter 
directly. The Matlab code includes multipath effect for various cut-off frequencies. 
The scope of this thesis does not cover the multipath effect, so multipath effect 
has been removed from the code and all the alphas are considered as one. 

Ezzat’s solution does not give an error ellipsoid for position estimate error, 
nor does it yield emitter position-covariance as a function of TDOA measurement 
error. 

The simulation of this thesis is based on Ezzat’s simulation. In addition to 
his simulation, position estimate error is based on multiple simulation runs to get 
position estimate covariance. 

The code of simulation is modified as in Appendix B to decide the number 
of experiments for minimum estimation error. The total standard deviation of the 
estimated position for multiple experiments is computed in Equation (3.77). The 
result is as in Figure 25. 
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Figure 25. The Effect of the Number of Experiments on Accuracy 


It can be seen from Figure 25 that the total standard deviation reaches a 
stable level around 500 experiments. 

The covariance matrix is based on the error resulting from the noise. 

Ezzat’s technique does not show the effect of the noise. The effect of the 
noise to TDOA has been added to this technique and simulation as in 


TDOA = t + e N 


(3.60) 


where 

t is the actual time difference of arrival 
e N is the error caused by the noise (Noise Error) 
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e w is a random number from a normal distribution with a 0 mean and <j t 

standard deviation. Standard deviation is explained in detail in Section C Part 2 
and the values used for the simulation are taken from that part. Because 
bandwidth, SNR and integration time changes for every receiver pair, noise error 
is different for each receiver pair, although the TDOA standard deviation is 
assumed to be the same for each receiver pair. 

The simulation is done for 500 times for efficiency as seen in Figure 25. 
Each result of the experiments for emitter location error is collected in an error 
matrix with a size of 500*3. 500 represents the number of experiments and 3 
represents the location estimation error for each axis. This error matrix is used to 
compute the covariance matrix by using Matlab. 

Covariance matrix can be found as 


cov(X) = E[(X - E[X])(X - E[X]) T ] (3.61) 


If Equation (3.61) is applied to multiple experiments, the covariance matrix 
consists of variances and correlation coefficients for every axis as seen in 
Equation (3.62). 


'uu 


a x Pxy°x°y Pxz < 7 x (J z 

Pxy^x^y ®y Pyz®y ^z 

Pxz°x°z PyzVy^z G z 


(3.62) 


The covariance matrix shows the relationship between every axis (x, y, 
and z). Because the covariance matrix shows the relationship between all axis, 
we convert to a covariance matrix which shows the only the variance for each 
single axis as in 
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(3.63) 


The covariance matrix in Equation (3.63) can be interpreted as a 
confidence ellipsoid. Standard deviation in every axis of Equation (3.63) 
represents half of the length of each axis of confidence ellipsoid. 2 <j x represents 

the length of the x-axis, 2 a y represents the length of the y-axis, and 2cr z 
represents the length of the y-axis of the ellipsoid. 

Because the covariance matrix in Equation (3.62) shows the noise 
correlation between every axis, it has to be converted to a diagonalized 
covariance matrix as in Equation (3.63) in order to get the relationships just for 
axes with themselves, which can be interpreted as standard deviations along 
three orthogonal position axes. 

Matlab command is used to get the eigenvalues of the covariance matrix. 
The eigenvalues matrix is the rotation matrix, which can be applied to covariance 
matrix in Equation (3.62) to get the diagonalized covariance matrix in Equation 

(3.63) . The following computation is done to apply the rotation matrix to the 
covariance matrix (3.62) in order to get the diagonalized covariance matrix in 

(3.63) 


DiagonalCovMat = (Eigenvectors)' * (CovMat)' 1 * Eigenvectors (3.64) 

The rotation matrix is explained in Slabaugh and Shoemake’s paper 
(Slabaugh, 1999 & Shoemake, 1985). The rotation has been done according to 
the following order; rotate x axis, rotate y axis and rotate z axis. 
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R=R z (^)R y (#)R x (^) 


( 3 . 65 ) 


where 


R z (^) is the rotation in z axis with the angle of (f> 

R y (#) is the rotation in y axis with the angle of # 

R x (^) is the rotation in y axis with the angle of y 


Rotation for z axis is defined as follows 


R Z W= 


cos^ 

sin^ 


0 


-sin^ 0 
cos^ 0 
0 1 


Rotation for y axis is defined as follows 


cos# 


R y (*)= 


0 

sin# 


0 

1 

0 


sin# 

0 

cos# 


Rotation for x axis is defined as follows 


R*(vO= 


i 

o 

0 


0 

cos \ff 
sin^ 


0 

-siny/ 
COS If/ 


(3.66) 


(3.67) 


(3.68) 
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As a result rotation matrix is 


R u /-? 12 /-? 13 

^21 ^22 ^23 


R 


31 


O P 

' '32 ' *33 


(3.69) 


and in terms of rotation angle 


R = 


cos 6 cos <j> 
costfsin^ 
-sin 0 


sin^sin0cos^-cos^sin0 
sin y/ sin 9 sin cj) + cos y/ cos <f> 
sin^cos^ 


cos^sin^cos^ + sin^sin^ 
cos y/ sin 9 sin (f>- sin y/ cos <j> 
cos ^ cos <9 


(3.70) 


where the rotation angles are 


<9 =-sin ~\R 3 ,) 


(3.71) 


f 

y/ = atan2 

v 


p R \ 
r *32 r *33 

cos#’cos 6 , 


(3.72) 


<t> = atan2 


» 2 , i q, 1 

v COS#’COS# y 


(3.73) 


The rotation matrix has to be orthogonal; in other words, R*R T =1, where 

I is the identity matrix or R 1 =R T , where R ’ 1 is inverse of rotation matrix and R T 
is the transpose of the rotation matrix. The easiest way to find the rotation matrix 
is to find the eigenvectors of the covariance matrix. A matrix formed from the 
combination of the eigenvectors of the covariance matrix is the rotation matrix. 

As explained earlier, the result in Equation (3.63) gives the variances for 

each axis diagonally. The square root of these variances gives the standard 
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deviation that is essential to check the accuracy of the system. An ellipsoid 
represents where the estimated emitter location is; the estimated location of the 
emitter might be any point inside that ellipsoid. These variances form this 
ellipsoid. Using the square root of each variance, standard deviations represent 
the length of the half of each axis. This ellipsoid represents the volume of space 
which contains the emitter location with probability 0.25. 

A 25 percent probability represents one standard deviation. We need to 
multiply these standard deviations to change the scale of the ellipsoid in order to 
add different probabilities. Ezzat’s algorithm and simulation do not include any 
probability of the estimated emitter location. 

The probability is added to the simulation as the length of the each semi¬ 
axis of the ellipsoid as explained in Torrieri’s paper as (Torrieri, 1984) 


Ix-axis ~ 4 K(J x 

= yfc J < 3 - 74 > 

4-„s= 


The probability of the actual location of the emitter lies inside that ellipsoid 
is given by (Torrieri, 1984) 


P e (rc) = erf(4K/2)-(42K/n)exp(-Kl2) (3.75) 


where 
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(3.76) 


erf(x) 


2 



Jexp {-t 2 )dt 


0 


If the reader is more interested in the statistical theory, he/she can read 
Torrieri’s paper (Torrieri, 1984). 


k- values have been calculated for various probabilities and shown in Table 


1 . 


Probability 

K 

99 % 

13.89702713 

95% 

8.934259824 

90 % 

6.922544695 

85% 

5.766298681 

80 % 

4.949459443 

75% 

4.314831079 

50 % 

2.298290951 

45% 

2.008418735 

25% 

1.011546612 


Table 1. Probability of Actual Location of Emitter to Lie in Ellipsoid and k 

Values 


The k values become larger when the probability increases, resulting in a 
larger ellipsoid as seen in Figure 26. 



Figure 26. k Values 
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C. SIMULATION RESULTS 

Since the standard deviation gives a good understanding of the accuracy, 
the simulation is run to analyze the accuracy with the standard deviation for 
every axis and the total standard deviation. 

Total standard deviation is given by 


a total = ^ a x 2 + G y + G z (3-77) 

where 

o x is the standard deviation for x-axis 

cr y is the standard deviation for y-axis 

o z is the standard deviation for z-axis 

The simulation gives the following parameters as output 

• Standard deviation for each axis as follows (in meters) 

Std_Dev_Data = 

29.2089 (x axis) 

30.0181 (y axis) 

80.2068 (z axis) 

• Total standard deviation as follows (in meters) 

T ot_Std_Dev = 

90.4842 

• Plot of the locations of the receivers are shown in Figure 27 and an 
exaggerated confidence ellipsoid is added in Figure 28. 
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Figure 27. TDOA Geometry in 3D 


The exaggerated confidence ellipsoid as in Figure 28 is the diagonalized 
confidence ellipsoid with axes exaggerated by a factor of 1000 for visibility and 
centered on the estimated position. 
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Figure 28. TDOA Geometry and Exaggerated Confidence Ellipsoid 

The simulation is run for the following analysis 

• The effect of the distance on the accuracy, 

• The effect of the noise on the accuracy, 

• Comparison between moving receivers and the stationary 
receivers, 

• The optimum distribution of the receivers for better accuracy. 

1. The Effect of the Distance on the Accuracy 

The developed simulation is in Appendix C. 

The effect of the distance on the accuracy is analyzed in different ways; 
distance change in one axis, distance changes in two axes, and distance 
changes in all three axes. Figure 29 shows the effect of the distance in one axis, 
Figure 30 shows the effect of the distance in two axes and Figure 31 shows the 
effect of the distances in three axes. 
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Figure 29. Effect of Distance to Accuracy for x-axis 
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Figure 30. Effect of Distance to Accuracy for two Axes (x and y axis) 
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Figure 31. Effect of Distance to Accuracy for Three Axes 


The distance affect the total standard deviation linearly since the distance 
change happens in one or two axes since these kind of changes result in a 
change in the shape of the geometry of the receiver location. The distance 
change in two axes affect the total standard deviation more rapidly compared to 
the distance change in one axis. 

On the other hand, total standard deviation is not effected linearly from 
distance change in three axes because changes in three axes result in changing 
the scale of the cluster, nothing more. The average total standard deviation in 
three axes distance change is around 300 meters for this configuration of random 
geometry. 

2. The Effect of the Noise on the Accuracy 

Noise affects the system and the simulation as 
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TDOA -t+e N (3.78) 

where t is the time of arrival and e w is the noise error. 

The noise error changes for each receiver pair. For this simulation, the 
noise distribution is considered as a normal distribution with a standard deviation 
<j x and zero mean for noise error. 

The developed simulation for stationary receivers is in Appendix D. 

The simulation is run for various standard deviations to understand the 
effect of the noise on the system. Figure 32 shows the effect of the noise on the 
accuracy. 



Figure 32. The Effect of Noise on Accuracy 


From Figure 32, it can be seen that the accuracy in terms of total standard 
deviation changes linearly with respect to noise standard deviation. This is a 
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good method of increasing the Signal to Noise Ratio (SNR) to keep the accuracy 
high. 

This is consistent with the theoretical result that the diagonalized 
variances of the estimated position are proportional to the measured TDOA 
variances. 



Figure 33. The Effect of Noise to Accuracy with Moving UAV 

From Figure 33, noise effect with moving UAV can affect the accuracy just 
like the way it affects without the UAV, but total standard deviation for estimated 
emitter location with UAV is smaller than the one without the UAV. 

3. Comparison Between Moving Receivers to the Stationary Receivers 

Ground forces can use UAVs instead of stationary receivers even though 
the use of stationary receivers is more probable. The simulation is developed to 
show the effect of moving UAVs to the accuracy in terms of total standard 
deviation. 


78 







The developed simulation is in Appendix F for both one and two flying 

UAVs. 

The simulation is run once for five random stationary receivers, four 
random stationary receivers plus a moving UAV and three random stationary 
receivers plus two moving UAVs. The simulation is run for random stationary 
receivers first, then the same procedure is followed for four random stationary 
receivers and a moving UAV. The same procedure is followed for the third 
scenario where there are three random stationary receivers plus two moving 
UAVs. For each set of receiver positions, the emitter is placed at (0,0,0) and 500 
TDOAs with random noise are generated. 

The results are analyzed in terms of the standard deviation for each axis 
and total standard deviation, which is a combination of standard deviations of ail 
axes. 

Figure 34 shows the random locations of five stationary receivers which 
try to locate the emitter. 



Figure 34. Random Locations of Stationary Receivers 
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Figure 35 shows the magnified confidence ellipsoid for five stationary 


receivers. 



Figure 35. Magnified Confidence Ellipsoid for Random Stationary Receivers 

Standard deviations for each axis for five random stationary receivers are: 
(t x = 30.3517 m 
o y = 32.7943 m 
<r z = 77.5429 m 

Total standard deviation for five random stationary receivers is: 

(7 T = 89.4964 m 

Figure 36 shows the random locations of four receivers and a moving UAV 
which try to locate the emitter. 
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The UAV flies over the operation area at a constant altitude for that 
scenario. 20 points of the route of the UAV are used for experiments. 500 
experiments are done at every point of the UAV. Total standard deviation is 
calculated for every point of the UAV. The results are collected in a matrix and 
plotted as in Figure 37. 


x 10 


TDOA Geometry 



Figure 36. Random Locations of Four Stationary Receivers and One Moving UAV 
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Figure 37. Total Standard Deviation Change for Four Random Stationary Receivers 

and a Moving UAV 

The confidence ellipsoid for each experiment is not plotted but the numeric 
values of standard deviation for each axes can be found in the proposed 
simulation code. 

Figure 38 shows the random locations of three stationary receivers and 
two moving UAV, which try to locate the emitter. 

The UAVs fly over the operation area at a constant altitude for that 
scenario. 20 points of the routes of each UAV are used for experiments. 500 
experiments are done at every point pair of the UAVs. Total standard deviation is 
calculated for every point pair of the UAVs. The results are collected in a matrix 
and plotted as in Figure 39. 
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Figure 38. Random Locations of Three Stationary Receivers and Two Moving UAVs 



Figure 39. Total Standard Deviation Change for Three Random Stationary Receivers 

and two Moving UAVs 
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The confidence ellipsoid for each experiment is not plotted but the numeric 
values of standard deviation for each axes can be found in the proposed code. 

It is clear from these datum that four random stationary receivers and a 
flying UAV has the smaller standard deviation for each axis, and they also have 
smaller total standard deviation compared to the other scenarios. 

4. The Optimum Distribution of the Receivers for Better Accuracy. 

Four receivers can be positioned around the emitter using different 
geometric patterns, such as straight, trapezoidal, parallelogram, inverted triangle, 
Y shaped, lozenge, square, and rectangle as explained in Yan-Ping’s article. 3D 
curves of position error are showed in Figure 40 and Figure 41. 




Figure 40. 3D Curves of Position Error for Straight (a), Trapezoidal (b), Parallelogram 

(c) and Rectangle (d) 


84 



















' e) Square 


(t) Y-E-harp' 


Figure 41. 3D Curves of Position Error for Lozenge (e), Inverted Triangle (f), Square 

(g), Y-Shaped (h) 


Coverage areas are shown in Table 2. 


Scene 

Coverage value 

Straight 

10.2% 

Trapezoidal 

29.4% 

Parallelogram 

40.5% 

Rectangle 

40.9% 

Lozenge 

54.9% 

Inverted triangle 

58.6% 

Square 

61.6% 

Y-sharp 

88.6% 


Table 2. Coverage Values for Different Type of Receiver Distributions (From 

Yan-Ping et al, 2010) 
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It is clear that the Y-shaped receiver distribution is the best for coverage 
values of the desired reach area. 

The developed simulation is in Appendix F. 

The following position parameters are used in the simulation for eight 
different scenarios, 

Scenario 1: Receivers are around the emitter (without altitude difference) 
as in Figure 42. The units are in kilometers (km). 


XI = 2; 

Y1 =0; 

Z1 =1.97; 

X2 = 0.7; 

Y2 = 2; 

Z2 = 1.98; 

X3 =-2; 

Y3 = 1.5; 

Z3 = 1.96; 

X 

II 

1 

M 

Y4 = -1.5; 

Z4 = 2; 

X5 = 1; 

Y5 = -1.5; 

Z5 = 1.99; 


Geometry and Confidence Ellipsoid 

1500 

1000 
500 

□ 

-500 

-1000 
-1500 


Figure 42. TDOA Geometry for Scenario 1 
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Scenario 2: Receivers are on a straight line (without altitude difference) as 
in Figure 43. The units are in kilometers (km). 


XI = -2; 

Y1 

= -1 

; zi 

= 

0.99; 

X2 = -1; 

Y2 

= 0 

Z2 

= 

i; 

X3 = 0; 

Y3 

= 2 

Z3 

= 

0.98; 

X4= 1; 

Y4 

= 3 

Z4 

= 

0.97; 

X5= 1.5; 

Y5 

= 4 

Z5 

— 




Scenario 3: Receivers are around the emitter (with altitude difference) as 
in Figure 44 and Figure 45. The units are in kilometers (km). 


XI = 2; 

Y1 =0; 

ZI =C 

X 

ro 

ii 

o 

cvf 

ii 

CM 

X 

Z2 = 1 

X3 =-2; 

Y3 = 1.5; 

Z3 = - 

X 

Ji¬ 

ll 

1 

ro 

Y4 = -1.5; 

Z4 = 2 

X5= 1; 

Y5 = -1.5; 

Z5 = - 
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TDOA Geometry and Confidence Ellipsoid 



Figure 44. TDOA Geometry for Scenario 3 (Top View) 


TDOA Geometry and Confidence Ellipsoid 



Figure 45. TDOA Geometry for Scenario 3 (Side View) 
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Scenario 4: Receivers are on a straight line (with altitude difference) as in 
Figure 46. The units are in kilometers (km). 


X 

II 

1 

M 

i 

ii 

X 

Z1 =0 

1 

II 

CM 

X 

o 

ii 

CM 

X 

Z2 = 1 

X3 = 0; 

Y3 = 2; 

Z3 = 2 

ii 

X 

Y4 = 3; 

Z4 = 3 

X5= 1.5; 

Y5 = 4; 

Z5 = 4 


TDOA Geometry and Confidence Ellipsoid 



Figure 46. TDOA Geometry for Scenario 4 

Scenario 5: Receivers are in a trapezoidal formation as in Figure 47. The 
units are in kilometers (km). 
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Scenario 6: Receivers are in a parallelogram formation as in Figure 48. 
The units are in kilometers (km). 
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Figure 48. TDOA Geometry for Scenario 6 

Scenario 7: Receivers are in a lozenge formation as in Figure 49. The 
units are in kilometers (km). 


XI = 

8; 

Y1 = 15; 

Z1 = 

-0.12 

X2 = 

-13; 

Y2 = 12.5; 

Z2 = 

-0.13 

X3 = 

0; 

Y3 = 18; 

Z3 = 

-0.14 

X4 = 

13; 

Y4 = 12.5; 

Z4 = 

-0.11 

X5 = 

0; 

-< 

cn 

ii 

00 

Z5 = 

-0.1; 
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Figure 49. TDOA Geometry for Scenario 7 


Scenario 8: Receivers are in an inverted triangle formation as in Figure 50. 


in kilometers 

(km). 


XI =2; 

Y1 =0; 

Z1 =0; 

X2 = 4; 

Y2 = 2; 

Z2 = -0 

X3 = 6; 

Y3 = -2; 

Z3 = -0 

X4 = 8; 

Y4 = 4; 

Z4 = -0 

X5 = 10; 

Y5 = -4; 

Z5 = -0 
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Table 3 shows the total standard deviations for all eight scenarios. 


Scenario 

Total Standard 

# 

Deviation 

1 

3.9003 km 

2 

1.2128 km 

3 

67.2223 m 

4 

751.0645 m 

5 

9.5230 km 

6 

1.4135 km 

7 

1.2795 km 

8 

3.7662 km 


Table 3. Estimated Emitter Location and Total Standard Deviation for Position 

Scenarios 
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It is clear from Table 3 that scenario 3 has the best two total standard 
deviations. This type of receiver positioning should be used to get the best 
accuracy if the system is using the proposed close-form geolocation technique. 
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V. CONCLUSION AND RECOMMENDATIONS 


A. CONCLUSION 

Electronic warfare plays a dominant role in today’s technological military 
world. Electronic warfare support is one of the core elements of EW, to include 
direction finding and geolocation. 

Passive geolocation of an emitter, requiring no emission of energy by the 
receiving and locating devices, plays an important role in both military and civilian 
applications. Passive geolocation can locate targets by measuring their 
electromagnetic emissions. 

The greatest internal problem for the Turkish military at present is the 
terrorist groups operating within the southeast part of Turkey. Military forces have 
been conducting operations against them for nearly 30 years. The operations 
occur in two phases, find and destroy. Finding the location of the terrorists in the 
mountainous areas of Turkey is crucial. Civilian and military intelligence has been 
used at the strategic level to locate them, but without operational and tactical 
location systems success is fleeting at best. 

One approach that uses multiple receivers is TDOA, which is the main 
concern of this thesis. For this technique, the time of arrival of the received signal 
at different receivers is measured; the measurements of the times of arrival are 
then used to find the time difference of arrival. These differences are used to 
determine the location of the emitter. 

Finding the time difference of arrival for radar signals is relatively easy, 
compared to CW signals. Because radar signals contain pulses, and these 
pulses are almost square (ideally), they have a defined beginning and an end. 
These beginnings and endings provide a time mark for time difference of arrival. 
These time marks are used for measuring the times of arrival; subtracting times 
of arrival results in a time difference of arrival. For a CW signal, there are no time 
marks. They do not have a specific beginning or end, which makes time 
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difference of arrival more difficult to calculate. A cross-correlation technique has 
to be used for CW signals to find the time difference of arrival. 

There is quite a bit of literature discussing many proposed methods to find 
the location of the emitter using the measured time difference of arrival. Many of 
the techniques require range data. The technique discussed in this thesis, 
Ezzat’s technique, does not. His technique is a closed-form geolocation 
technique. 

The areas where the operations take place are very mountainous. 
Because of this reason, it is necessary to use a 3D TDOA technique. 

This thesis provided a brief discussion on the fundamentals of TDOA and 
a closed-form usage of TDOA. A simulation model was developed in order to 
observe the performance of the proposed close-from technique in different 
scenarios. The effects of the distance, noise and distribution of the receivers to 
accuracy were analyzed. A comparison between a system of stationary receivers 
and a system of combinations of stationary and moving receivers was done to 
understand the best combination of receivers for best accuracy. For each 
combination of receivers, 500 simulations with added independent noise were 
run to acquire an estimate of the emitter position and to obtain diagonalized 
covariance values. 

The developed simulation was based on Ezzat’s closed-form technique. 
His technique was repeated many times to understand its robust nature. 

Distance affects accuracy in a linear manner, as explained in Chapter 4 
Section C. Being close to the emitter gives the best accuracy. The closer to the 
emitter, the better the accuracy. 

Noise also has a linear effect on accuracy, as explained in Chapter 4 
Section C. When the standard deviation of the noise increases, the standard 
deviation of the position estimate increases proportionally, that results a 
reduction in the accuracy. Since the noise considered in the simulation is 
background noise and cannot be removed from the system, the best way to 
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overcome this problem is to use a receiver with a better sensitivity so as to 
achieve a higher signal to noise ratio. 

The system should include at least four stationary and one flying UAV as 
the fifth receiver for better accuracy as explained in Chapter 4 Section C. It is 
clear from the simulation that having a flying UAV with stationary ground 
receivers has a significant effect on the accuracy. One other advantage of having 
a UAV is not losing the line of sight with the signal of interest. It is difficult to have 
line of sight at all times, especially in high clutter areas that the Turkish military 
forces operate. UAVs can solve this problem. They are also faster and more 
mobile than ground receivers, which make them preferable for better coverage. 
UAVs also offer the advantage of adding the FDOA technique to the calculation 
of the emitter/enemy location, a discussion that was beyond the interest area of 
this thesis, but is relevant to the objective of locating emitters. 

The two best distributions of receivers are shown in Scenarios 3 and 4, as 
explained in Chapter 4 Section C. One solution is to distribute the receivers 
around the emitter with altitude differences as in Figures 51 and 52. Figure 51 
shows the geometry from the side view and Figure 52 shows the geometry from 
the top view. The other distribution is to have the receivers on a baseline with 
altitude differences as in Figure 53. The operational forces should try to locate 
the receivers as explained in one of these scenarios for best accuracy. 
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TDOA Geometry and Confidence Ellipsoid 



Figure 51. TDOA Geometry of the Receivers around the Emitter (Top View) 


TDOA Geometry and Confidence Ellipsoid 



Figure 52. TDOA Geometry of the Receivers around the Emitter (Side View) 
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TDOA Geometry and Confidence Ellipsoid 
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Figure 53. TDOA Geometry of the Receivers on a Straight Line with Altitude 

Difference 



The main disadvantage of the proposed closed-form geolocation system is 
a requirement to use one more receiver than with some other TDOA techniques, 
for which it is only necessary to have four TDOA receivers. 

The simulation enables the user to analyze the performance of the closed- 
form geolocation technique under desired conditions. It must be noted that the 
values obtained from the simulation may differ from real environment values, 
since many assumptions were made and certain values were kept constant in 
order to simplify the optimization process. 

From the Turkish military perspective, since the developed model 
examines specific conditions, the model can be upgraded to an advanced level 
according to the needs and capabilities of the Turkish military and can be an 
example model for improved future studies. 
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B. RECOMMENDATIONS 


EW will maintain a critical place on the battlefield. Almost all nations are 
spending money for EW research. Every nation tries to take part in the 
intelligence arena. All these facts increase the importance of operational tactics 
development for these technologies. 

It is obvious that closed-form geolocation has certain advantages in 
specified areas. The research conducted under this thesis tried to examine these 
features and created a model in order to simulate the environment. The 
simulation model can be a starting point and can be improved in several ways for 
future studies. 

Contrasting the proposed closed-form geolocation technique with the 
other forms of TDOA geolocation that use repeated measurements involving at 
least one moving receiver might be a good research subject to determine which 
TDOA geolocation technique to use for better accuracy. 

Instead of making assumptions, using fixed values and omitting certain 
facts in order to reduce the complexity, allowing for more detailed studies. 
Multipath effects can be added for better understanding. Allowing for multipath 
may well produce better results in the presence of mountainous clutter. The 
academic studies that have been used as references in this thesis are the 
primary resources that can lead to future research and improvements in TDOA 
geolocation. 

Future studies may try to answer how to decrease the number of receivers 
while maintaining accuracy. Because synchronization of the receivers plays a 
very important role in multiple receiver systems, as in this thesis, future studies 
may focus on how develop better synchronization in terms of system design and 
used techniques. 

A new design of the overall system may be another area of research, 
where the TDOA system and the receivers are connected to a central control 
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point over the Internet or over a tactical network. This may help decision makers 
to develop better battlefield awareness and more tactical options. 
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APPENDIX A 


MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 

TECHNIQUE 


function TDOA 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 

% Equations, % 2006 

% Latest revision: August 15, 2012 

%clear the screen 
clc 

%Variables 

%number of experiments for accuracy 
n = 500; 

%number of the positions of the UAV 
m = 1 ; 


k=8.934259824; %k for 95% of probability 


o, 

o 

k=6.922544695; 

%k 

for 

90% 

of 

probability 

o, 

o 

k=5.766298681; 

%k 

for 

85% 

of 

probability 

o, 

o 

k=4.949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

k=4.314831079; 

%k 

for 

75% 

of 

probability 


%Coordinates of emitter 
E=[0 0 0 ] ; 

%Speed of Light(Propagation) 
c = 3 * 10 A 8 ; 

%Set All Alphas to 1 

alphal = 1; 

alpha2 = 1; 

alpha3 = 1; 

alpha4 = 1; 

alpha5 = 1; 

%Stationary Receivers 
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X3 

=-1000; 

Y3 = 

2000; 

Z3 = 

-2000; 

X4 

= -1000; 

Y4 = 

-1000; 

Z4 = 

1000; 

X5 

= -3000; 

Y5 = 

-1000; 

Z5 = 

-1000; 


for i = 1:n, 

%initial coordinate matrix 
Coor2 = [0 0 0]; 

%Coordinates of Receivers 


for j = 

: 1: m. 





XI 

= 1000*j; 

Yl = 

-3000*j; 

Zl = 

1000; %Moving 

UAVs 

X2 

= 2000*j; 

Y2 = 

-2000*j; 

Z2 = 

-4000*j; 


%UAV Coordinate MAtrix 
XlCoor(1,j) = XI; 

YlCoor (1,j) = Yl; 

ZlCoor (1,j) = Zl; 

X2Coor(1,j) = X2; 

Y2Coor(1,j) = Y2; 

Z2Coor (1,j) = Z2; 

%Calculate the Distance From emitter to the 
Receiver for TOA 

ECl=sqrt(abs( (E(1)-XI) A 2+(E(2)-Yl) A 2+ (E (3)-Zl) A 2)) ; 
EC2=sqrt(abs( (E(1)-X2) A 2+(E(2)-Y2) A 2+(E(3)-Z2) A 2)) ; 
EC3=sqrt(abs( (E(1)-X3) A 2+(E(2)-Y3) A 2+ (E (3)-Z3) A 2)) ; 
EC4=sqrt(abs((E(1)-X4) A 2+(E(2)-Y4) A 2+(E(3)-Z4) A 2)); 
EC5=sqrt (abs( (E(1)-X5) A 2+(E (2)-Y5) A 2+(E(3)-Z5) A 2)) ; 

%Generate random noise for every channel related to 
Standat Deviation 

Noise Err = random(' Normal 0,5e-8, 1,5) ; 


%Calculate 

the times 

of 

arrival 


tl 

= EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err(1); 

t2 

= EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

~Err (2) ; 

t3 

= EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err(3) 

14 

= EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err (4) 

t5 

= EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err (5); 


%from now on calculate the emitter location 
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%Calculate the coefficients in Eqs.(17). Note that 
alpha (assumed),not alpha actual is used. 

all = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t3-tl) *(X3/alpha3 A 2 - XI/alphal A 2); 

al2 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t3-tl) * (Y3/alpha3 A 2 - Y1/alphal A 2); 

al3 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t3-tl) * (Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 

((X3 A 2+Y3 A 2+Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (t3-t2); 

a21 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/ (14-t1) *(X4/alpha4 A 2 - XI/alphal A 2); 

a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/ (14-t1) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t4-tl) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (14-12) ; 

a31 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t5-tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t5-tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t5-tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

((X5 A 2+Y5 A 2+Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (15-12) ; 

%Calculate the matrix A on the left-hand side. 

A = double([[all al2 al3]; [a21 a22 a23]; [a31 a32 

a 3 3 ] ] ) ; 


%Calculate the vector B on the right-hand side. 
B = double([bl; b2; b3]); 


%Solve the simultaneous equations AX = B. 
Coorl = (A'*A)\A'*B; 

Coor2 =+ Coorl; 
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end 


Coor = Coor2/m; 

%Error by every Coordinate 

Err = [abs(Coor (1)-E(1)),abs (Coor (2)-E (2)), 
abs (Coor(3)—E(3))]; 

%Error Matrix 
Mat(i,:) = Err; 

%Total Error 

Total_Err = sqrt(Err(1) A 2 + Err(2) A 2 + Err(3) A 2); 


end 

%calculate Cov.Matrix 
Cov_Mat = cov(Mat); 

%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues ] = eig(Cov_M a t); 

%Find angles (from "Computing Euler angles from a rotation 
matrix") 

Beta = -asind (Eigenvectors (3, 1) ) ; 

Gamma = 

(atan2 (Eigen_V ec tors(2,1)/cosd(Beta),Eigen_V ec tors(1, 1) / cos 
d(Beta)))*(180/pi); 

Alpha = 

(atan2(Eigen_V ec tors(3,2)/cosd(Beta),Eigen_V ec tors(3,3)/cos 
d(Beta)))*(180/pi); 

%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*k) 

%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 

%Draw the error ellipsoid from the diagonalized covariance 

matrix 

hold on 

grid on 
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%ellipsoid 
[x, y, z] = 

ellipsoid(E(1),E(2),E(3),Std_Dev_Data(1),Std_Dev_Data(2),St 
d_Dev_Data(3),20); 

%color of the ellisoid 
colormap copper; 

%ploy ellipsoid 
hMesh = mesh(x,y,z); 

%rotate ellipsoid 
rotate(hMesh,[1 0 0],Alpha); 
rotate(hMesh,[0 1 0],Beta); 
rotate(hMesh,[0 0 1],Gamma); 

%equal axis 
axis equal 

%# Change the camera viewpoint 
view([-36 18]); 

%label axises and title 
xlabel( 'X-Axis' ) 
ylabel( 'Y-Axis' ) 
zlabel( 'Z-Axis' ) 

title( 'TDOA Geometry and Confidence Ellipsoid') 

%plot stationary receivers 

plot3 (X3,Y3,Z3, 'o' ); 
plot3 (X4,Y4,Z4, 'o' ); 
plot3 (X5,Y5,Z5, 'o' ); 

%plot moving receiver 

plot3(XlCoor,YlCoor,ZlCoor, 'ro' ); 

plot3(X2Coor,Y2Coor,Z2Coor, 'go' ); 

%plot emitter 

plot3 (E (1),E (2),E (3), 'ms' ); 
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APPENDIX B 


MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 
TECHNIQUE TO DECIDE THE NUMBER OF EXPERIMENTS FOR OPTIMUM 

ACCURACY. 


function TDOA_Num_Experiment_vs_Accurracy 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 

% Equations, % 2006 

% Latest revision: August 15, 2012 

%clear the screen 
clc 

%Variables 

%number of experiments for accuracy 
n = 700; 

%number of the positions of the UAV 
m = 1 ; 


kk=8.934259824; 

% k=6.922544695; 
% k=5.766298681; 
% k=4.949459443; 
% k=4.314831079; 


%k for 95% of probability 


%k 

for 

90% 

of 

probability 

%k 

for 

85% 

of 

probability 

%k 

for 

80% 

of 

probability 

%k 

for 

75% 

of 

probability 


%Coordinates of emitter 
E=[0 0 0 ] ; 


%Speed of Light(Propagation) 
c = 3 * 10 A 8 ; 


%Set All Alphas to 1 

alphal = 1; 

alpha2 = 1; 

alpha3 = 1; 

alpha4 = 1; 

alpha5 = 1; 


%Stationary Receivers 
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XI = 

1000; 

Y1 

X2 = 

1000; 

Y2 

X3 = 

-1000; 

Y3 

X4 = 

-1000; 

Y4 

X5 = 

-3000; 

Y5 

for 

k = 2 : n. 



for i = 2:k. 


1000; 

Z1 

= 1000; 

2000; 

Z2 

= 3000; 

2000; 

Z3 

= -2000; 

-1000; 

Z4 

= 1000; 

-1000; 

Z5 

= -1000; 


%Calculate the Distance From emitter to the 
Receiver for TOA 


EC1 = 

=sqrt 

(abs 

(E 

(1) 

-XI) 

A 2 + 

(E 

(2) 

-Yl) 

A 2 + 

(E 

(3) 

-Zl) 

A 2) ) ; 

EC2 = 

=sqrt 

(abs 

(E 

(1) 

-X2) 

A 2 + 

(E 

(2) 

-Y2) 

A 2 + 

(E 

(3) 

-Z2) 

A 2) ) ; 

EC3= 

=sqrt 

(abs 

(E 

(1) 

-X3) 

A 2 + 

(E 

(2) 

-Y3) 

A 2 + 

(E 

(3) 

-Z3) 

A 2) ) ; 

EC4 = 

=sqrt 

(abs 

(E 

(1) 

-X4) 

A 2 + 

(E 

(2) 

-Y4) 

A 2 + 

(E 

(3) 

-Z4) 

A 2) ) ; 

EC5= 

=sqrt 

(abs 

(E 

(1) 

-X5) 

A 2 + 

(E 

(2) 

-Y5) 

A 2 + 

(E 

(3) 

-Z5) 

A 2) ) ; 


%Generate random noise for every channel related to 
Standat Deviation 

Noise Err = random(' Normal 0,5e-8,1,5); 


%Calculate the times of arrival. 


tl 

= EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err(1) 

r 

t2 

= EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err (2) 

r 

t3 

= EC3 

/ 

(alpha3 


c) 

+ 

Noise 

"Err (3) 

r 

t4 

= EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err(4) 

r 

t5 

= EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err (5) 

r 


%from now on calculate the emitter location 

%Calculate the coefficients in Eqs.(17). Note that 
alpha (assumed),not alpha actual is used. 

all = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t3-tl) *(X3/alpha3 A 2 - XI/alphal A 2); 

al2 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t3-tl) *(Y3/alpha3 A 2 - Y1/alphal A 2); 

al3 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t3-tl) *(Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 

( (X3 A 2+Y3 A 2 + Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2 + Z1 A 2)/alphal A 2) + 
c A 2 * (t3-t2); 

a21 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/ (14-t1) *(X4/alpha4 A 2 - XI/alphal A 2); 
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a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/ (14-t1) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/ (14-t1) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (t4-t2); 

a31 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t5-tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t5-tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t5-tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

((X5 A 2+Y5 A 2+Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (t5-t2); 

%Calculate the matrix A on the left-hand side. 

A = double([[all al2 al3]; [a21 a22 a23]; [a31 a32 

a 3 3 ] ] ) ; 

%Calculate the vector B on the right-hand side. 

B = double([bl; b2; b3]); 

%Solve the simultaneous equations AX = B. 

Coor = (A'*A)\A'*B; 

%Error by every Coordinate 

Err = [abs(Coor (1)-E (1)),abs (Coor (2)-E (2)), 
abs (Coor (3)-E (3))]; 

%Error Matrix 
Mat (i, :) = Err; 


end 

%calculate Cov.Matrix 
Cov_Mat = cov(Mat); 

%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues ] = eig (Cov_M a t) ; 

%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*kk) ; 
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%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 


Tot Std Devi(k,:) = Tot Std Dev; 


end 

N = linspace (1,n,n); 

%plot the result 
plot (N,Tot_Std_Devl) 

%label axises and title 

xlabel (' Number of Experiments') 

ylabel (' Total Standard Deviation (m) ') 

title (' Effect of Number of Experiment to the Accuracy') 
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APPENDIX C 


MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 
TECHNIQUE TO SEE THE EFFECTS OF DISTANCE CHANGES ON 

ACCURACY. 

function TDOA_Distance_to_Accuracy 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 
% Equations, % 2006 
% Latest revision: August 15, 2012 

%clear the screen 
clc 

%Variables 


%number of experiments for accuracy 
n = 500; 


kk=8.934259824; %k for 95% of probability 


k=6.922544695; %k for 90% 
k=5.766298681; %k for 85% 
k=4.949459443; %k for 80% 
k=4.314831079; %k for 75% 


of probability 
of probability 
of probability 
of probability 


%Coordinates of emitter 
E=[0 0 0 ] ; 

%Speed of Light(Propagation) 
c = 3 * 10 A 8 ; 

%Set All Alphas to 1 

alphal = 1; 

alpha2 = 1; 

alpha3 = 1; 

alpha4 = 1; 

alpha5 = 1; 

%distance multipication constant 
k=10 0; 


%Distance Change loop 
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for 1 


1: k 


%Distance Change in three axis 


%Stationary Rece 

ivers 




XI 

= 1000*1; 

Y1 = 

1000*1; 

Z1 = 

= 1000*1; 

X2 

= 1000*1; 

Y2 = 

2000*1; 

Z2 = 

= 3000*1; 

X3 

=-1000*1; 

Y3 = 

2000*1; 

Z3 = 

= -2000*1 

X4 

= -1000*1; 

Y4 = 

-1000*1; 

Z4 = 

= 1000*1; 

X5 

= -3000*1; 

Y5 = 

-1000*1; 

Z5 = 

= -1000*1 


o, 

O 

Distance Change in two 

axis 




O, 

O 

%Stationary Receivers 





Q, 

O 

XI = 1000; 

Y1 = 

1000*1; 

Z1 

= 1000*1; 

O, 

O 

X2 = 1000; 

Y2 = 

2000*1; 

Z2 

= 3000*1; 

O, 

O 

X3 =-1000; 

Y3 = 

2000*1; 

Z3 

= -2000*1; 

o, 

o 

X4 = -1000; 

Y4 = 

-1000*1; 

Z4 

= 1000*1; 

Q, 

O 

O, 

X5 = -3000; 

Y5 = 

-1000*1; 

Z5 

= -1000*1; 

O 

O, 

O 

Distance Change in one 

axis 




O, 

O 

Stationary Receivers 





o, 

o 

XI = 1000; 

Y1 = 

1000; 

Z1 = 

1000*1; 

o, 

o 

X2 = 1000; 

Y2 = 

2000; 

Z2 = 

3000*1; 

Q. 

O 

X3 =-1000; 

Y3 = 

2000; 

Z3 = 

-2000*1; 

Q. 

O 

X4 = -1000; 

Y4 = 

-1000; 

Z4 = 

1000*1; 

O 

O 

X5 = -3000; 

Y5 = 

-1000; 

Z5 = 

-1000*1; 


for i = 2:n, 

%Calculate the Distance From emitter to the 
Receiver for TOA 


EC1 

=sqrt 

(abs ( 

(E (1 

) -xi) 

A 2 + 

(E 

(2) 

- Y1 

) A 2 + 

(E 

(3 

) -zi) 

A 2 ) 

) ; 

EC2 

=sqrt 

(abs ( 

(E (1 

) -X2) 

A 2 + 

(E 

(2) 

-Y2 

) A 2 + 

(E 

(3 

) -Z2) 

A 2 ) 

) ; 

EC3 

=sqrt 

(abs ( 

(E (1 

) -X3) 

A 2 + 

(E 

(2) 

-Y3 

) A 2 + 

(E 

(3 

) — Z 3) 

A 2 ) 

) ; 

EC4 

=sqrt 

(abs ( 

(E (1 

) -X4) 

A 2 + 

(E 

(2) 

-Y4 

) A 2 + 

(E 

(3 

) — Z 4 ) 

A 2 ) 

) ; 

EC5 

=sqrt 

(abs ( 

(E (1 

) -X5) 

A 2 + 

(E 

(2) 

-Y5 

) A 2 + 

(E 

(3 

) — Z 5) 

A 2 ) 

) ; 

%Generate random 

noise 

fo 

r 

eve 

ry 

channel 

related 

to 


Standat Deviation 

Noise Err = random(' Normal 0,5e-8,1,5); 


%Calculate 

the times 

of 

arrival 

. 



tl 

= EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err 

(1) 

r 

t2 

= EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err 

(2) 

r 

t3 

= EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err 

(3) 

r 

14 

= EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err 

(4) 

r 
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t5 = EC5 / (alpha5 * c) + Noise_Err(5); 


%from now on calculate the emitter location 


%Calculate the coefficients in Eqs.(17). Note that 
alpha (assumed),not alpha actual is used. 

all = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t3-tl) *(X3/alpha3 A 2 - XI/alphal A 2); 

al2 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t3-tl) *(Y3/alpha3 A 2 - Y1/alphal A 2); 

al3 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t3-tl) *(Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 

((X3 A 2+Y3 A 2+Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (t3-t2); 

a21 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/ (14-t1) *(X4/alpha4 A 2 - XI/alphal A 2); 

a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/ (14-t1) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/ (14-t1) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (14-12) ; 

a31 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t5-tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t5-tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t5-tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/ (t2-tl) * ((X2 A 2+Y2 A 2 + Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

( (X5 A 2+Y5 A 2 + Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2 + Z1 A 2)/alphal A 2) + 
c A 2 * (t5-t2) ; 

%Calculate the matrix A on the left-hand side. 

A = double([[all al2 al3]; [a21 a22 a23]; [a31 a32 

a 3 3 ] ] ) ; 


%Calculate the vector B on the right-hand side. 
B = double([bl; b2; b3]); 
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%Solve the simultaneous equations AX = B. 
Coor = (A'*A)\A'*B; 

%Error by every Coordinate 

Err = [abs(Coor (1)-E (1)),abs (Coor (2)-E (2)), 
abs (Coor (3)-E (3))]; 

%Error Matrix 

Mat (i, :) = Err; 


end 

%calculate Cov.Matrix 
Cov_Mat = cov(Mat); 

%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues ] = eig (Cov_M a t) ; 

%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*kk); 

%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 


Tot Std Devl(l,:) = Tot Std Dev; 


end 

N = linspace(1,k,k); 

%plot the result 
plot (N,Tot_Std_Devl) 

%label axises and title 

xlabel (' Distance Multiplier') 

ylabel (' Total Standard Deviation (m) ') 

title (' Effect of Distance to the Accuracy') 
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APPENDIX D 


MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 
TECHNIQUE WITH STATIONARY RECEIVERS TO SEE THE EFFECTS OF 

NOISE ON ACCURACY. 

function TDOA_Noise_vs_Accuracy_Stationary_Receivers 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 

% Equations, % 2006 

% Latest revision: August 15, 2012 

%clear the screen 
clc 

%Variables 

%number of experiments for accuracy 
n = 500; 

%standard deviation of noise chage constant 
1=2 5; 


% k values for probability 

kk=8.934259824; %k for 95% of probability 


o, 

o 

kk=6.922544695; 

%k 

for 

90% 

of 

probability 

o 

o 

kk=5.766298681; 

%k 

for 

85% 

of 

probability 

g, 

o 

kk=4.949459443; 

%k 

for 

80% 

of 

probability 

g, 

o 

kk=4.314831079; 

%k 

for 

75% 

of 

probability 


%Coordinates of emitter 
E=[0 0 0 ] ; 

%Speed of Light(Propagation) 
c = 3 * 10 A 8; 

%Set All Alphas to 1 

alphal = 1; 

alpha2 = 1; 

alpha3 = 1; 

alpha4 = 1; 

alpha5 = 1; 
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%Stationary 

Receivers 



XI 

= 1000; 

Yl 

= 1000; 

Zl = 

= 1000; 

X2 

= 1000; 

Y2 

= 2000; 

Z2 = 

= 3000; 

X3 

=-1000; 

Y3 

= 2000; 

Z3 = 

= -2000; 

X4 

= -1000; 

Y4 

= -1000; 

Z4 = 

= 1000; 

X5 

= -3000; 

Y5 

= -1000; 

Z5 = 

= -1000; 


%for loop for changing the standard deviation 
for k = 1:1, 

%for loop for experiment for every standard deviation 
for i = 1:n, 

%Calculate the Distance From emitter to the 
Receiver for TOA 

ECl=sqrt(abs( (E(1)-XI) A 2+(E(2)-Yl) A 2+(E(3)-Zl) A 2)) ; 
EC2=sqrt(abs( (E(1)-X2) A 2+(E (2)-Y2) A 2+ (E (3)-Z2) A 2)) ; 
EC3=sqrt(abs((E(1)-X3) A 2+(E(2)-Y3) A 2+(E(3)-Z3) A 2)); 
EC4=sqrt(abs((E(1)-X4) A 2+(E(2)-Y4) A 2+(E(3)-Z4) A 2)); 
EC5=sqrt (abs ( (E (1) -X5) A 2+ (E (2) -Y5) A 2+(E(3)-Z5) A 2)); 

%Generate random noise for every channel related to 
Standat Deviation 

Noise Err = random(' Normal 0,k*le-8, 1,5) ; 


%Calculate 

the times 

of 

arrival 

• 

tl 

= EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

_Err (1) 

t2 

= EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err (2) ; 

t3 

= EC3 

/ 

(alpha3 

-k 

c) 

+ 

Noise 

“Err (3) 

t4 

= EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err (4) 

t5 

= EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err (5); 


%from now on calculate the emitter location 

%Calculate the coefficients in Eqs.(17). Note that 
alpha (assumed),not alpha actual is used. 

all = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t3-tl) *(X3/alpha3 A 2 - XI/alphal A 2); 

al2 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t3-tl) * (Y3/alpha3 A 2 - Yl/alphal A 2); 

al3 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t3-tl) *(Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 
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((X3 A 2+Y3 A 2+Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (13-12); 

a21 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/ (14-t1) *(X4/alpha4 A 2 - XI/alphal A 2); 

a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/ (14-t1) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/ (14-t1) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (14-12) ; 

a31 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t5-tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t5-tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t5-tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

( (X5 A 2+Y5 A 2 + Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2 + Z1 A 2)/alphal A 2) + 
c A 2 * (t5-t2) ; 

%Calculate the matrix A on the left-hand side. 

A = double([[all al2 al3]; [a21 a22 a23]; [a31 a32 

a 3 3 ] ] ) ; 


%Calculate the vector B on the right-hand side. 
B = double([bl; b2; b3]); 

%Solve the simultaneous equations AX = B. 

Coor = (A'*A)\A'*B; 

%Error by every Coordinate 

Err = [abs(Coor(1)-E(1)),abs (Coor (2)-E (2)), 
abs (Coor (3)-E (3))]; 

%Error Matrix 
Mat(i,:) = Err; 

end %i 

%calculate Cov.Matrix 
Cov Mat = cov(Mat); 
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%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues] = eig (Cov_M a t) ; 

%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*kk); 

%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 


Tot Std Devi(k,:) = Tot Std Dev; 


end %k 


N = linspace (1,1,1); 


%plot the result 
plot (N,Tot_Std_Devl) 


%label axises 
xlabel( 'Noise 
ylabel( 'Total 
title ( 'Effect 


and title 
Standard 
Standard 
of Noise 


Deviation (sec) ') 
Deviation (m)' ) 
to the Accuracy') 


120 


APPENDIX E 


MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 
TECHNIQUE WITH ONE MOVING UAV TO SEE THE EFFECTS OF NOISE ON 

ACCURACY. 

function TDOA_Noise_vs_Accuracy_With_Flying_UAV 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 

% Equations, % 2006 

% Latest revision: August 15, 2012 

%clear the screen 
clc 

%Variables 

%number of experiments for accuracy 
n = 500; 

Mat = zeros(n,3); 

%number of the positions of the UAV 
m = 10; 

%standard deviation of noise chage constant 
1=2 5; 


% k values for probability 

kk=8.934259824; %k for 95% of probability 


o, 

o 

kk=6.922544695; 

%k 

for 

90% 

of 

probability 

o, 

o 

kk=5.766298681; 

%k 

for 

85% 

of 

probability 

o 

o 

kk=4.949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

kk=4.314831079; 

%k 

for 

75% 

of 

probability 


%Coordinates of emitter 
E=[0 0 0 ] ; 

%Speed of Light(Propagation) 
c = 3 * 10 A 8 ; 

%Set All Alphas to 1 
alphal = 1; 
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alpha2 = 1; 
alpha3 = 1; 
alpha4 = 1; 
alpha5 = 1; 


%Stationary 

Receivers 



X2 = 1000; 

Y2 

= 2000; 

Z2 = 

3000; 

X3 =-1000; 

Y3 

= 2000; 

Z3 = 

-2000; 

X4 = -1000; 

Y4 

= -1000; 

Z4 = 

1000; 

X5 = -3000; 

Y5 

= -1000; 

Z5 = 

-1000; 


%for loop for changing the standard deviation 
for k = 1:1, 


%for loop for experiment for every standard deviation 
for i = 1:n, 

%initial coordinate matrix 
Coor2 = [0 0 0]; 

%for loop for the positions of the UAV 
for j = l:m, 

XI = 1000*j; Y1 = 1000*j; Z1 = 1000; 

%Moving UAV 


Receiver 

Zl) A 2)); 
Z 2 ) A 2 ) ) ; 
Z3) A 2)); 
Z4 ) A 2)) ; 
Z5) A 2)) ; 


%Calculate the Distance From emitter to the 
for TOA 


EC1 = 

=sqrt 

(abs ( 

(E 

(1) 

-XI) 

A 2 + 

(E 

(2) 

-Yl) 

A 2 + 

(E 

(3)- 

EC2 = 

=sqrt 

(abs ( 

(E 

(1) 

-X2) 

A 2 + 

(E 

(2) 

-Y2) 

A 2 + 

(E 

(3)- 

EC3= 

=sqrt 

(abs ( 

(E 

(1) 

-X3) 

A 2 + 

(E 

(2) 

-Y3) 

A 2 + 

(E 

(3)- 

EC4 = 

=sqrt 

(abs ( 

(E 

(1) 

-X4) 

A 2 + 

(E 

(2) 

-Y4) 

A 2 + 

(E 

(3)- 

EC5= 

=sqrt 

(abs ( 

(E 

(1) 

-X5) 

A 2 + 

(E 

(2) 

-Y5) 

A 2 + 

(E 

(3)- 


%Generate random noise for every channel 
related to Standat Deviation 

Noise_Err = random (' Normal 0,k*le-8,1,5); 

%Calculate the times of arrival. 

tl = EC1 / (alphal * c) + Noise_Err(1); 

t2 = EC2 / (alpha2 * c) + Noise_Err(2); 
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t3 = EC3 / (alpha3 * c) + Noise_Err(3); 
t4 = EC4 / (alpha4 * c) + Noise_Err(4); 
t5 = EC5 / (alpha5 * c) + Noise_Err(5); 


%from now on calculate the emitter location 


%Calculate the coefficients in Eqs.(17). Note 
that alpha (assumed),not alpha actual is used. 

all = 2/ (t2-t1) * (X2/alpha2 A 2 - Xl/alphal A 2) 
2/(t3-tl) *(X3/alpha3 A 2 - XI/alphal A 2); 

al2 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) 
2/(t3-tl) * (Y3/alpha3 A 2 - Y1/alphal A 2); 

al3 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) 
2/(t3-tl) * (Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ( (X2 A 2+Y2 A 2 + Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 

((X3 A 2+Y3 A 2+Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (13-12); 

a21 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) 

2/ (14-t1) *(X4/alpha4 A 2 - XI/alphal A 2); 

a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) 

2/ (14-t1) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) 

2/ (14-t1) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/(t2-tl) * ( (X2 A 2+Y2 A 2 + Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (14-12) ; 

a31 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) 
2/(t5-tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) 
2/(t5-tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) 
2/(t5-tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

((X5 A 2+Y5 A 2+Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (t5-t2); 

%Calculate the matrix A on the left-hand side. 
A = double([[all al2 al3]; [a21 a22 a23]; [a31 

a32 a33]]); 


%Calculate the vector B on the right-hand side 
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B = double([bl; b2; b3]); 

%Solve the simultaneous equations AX = B. 
Coorl = (A'*A)\A'*B; 

Coor2 =+ Coorl; 
end % j 

Coor = Coor2/m; 

%Error by every Coordinate 

Err = [abs(Coor(1)-E(1)),abs (Coor (2)-E (2)), 
abs (Coor (3)-E (3))]; 

%Error Matrix 

Mat (i, :) = Err; 


end %i 

%calculate Cov.Matrix 
Cov_Mat = cov(Mat); 

%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues ] = eig (Cov_M & t) ; 


%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*kk); 

%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 


Tot_Std_Devl(k,:) = Tot_Std_Dev; 
end %k 

N = linspace (1,1,1); 

%plot the result 
plot (N,Tot_Std_Devl) 
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%label axises 
xlabel( 'Noise 
ylabel( 'Total 
title ( 'Effect 


and title 
Standard 
Standard 
of Noise 


Deviation (sec))') 
Deviation (m)' ) 
to the Accuracy') 
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APPENDIX F 


MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 
TECHNIQUE TO SEE THE EFFECTS OF HAVING A MOVING UAV ALONG 
WITH STAIONARY RECEIVERS ON ACCURACY. 

function TDOA_Flying_UAV 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 
% Equations, % 2006 
% Latest revision: August 26, 2012 

%clear the screen 
clc 

%Variables 

%number of experiments for accuracy 
n = 500; 

Mat = zeros(n,3); 

%number of the positions of the UAV 
m = 2 0; 


k=8.934259824; %k for 95% of probability 


o, 

o 

k=6.922544695; 

%k 

for 

90% 

of 

probability 

o 

o 

k=5.766298681; 

%k 

for 

85% 

of 

probability 

o, 

o 

k=4.949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

k=4.314831079; 

%k 

for 

75% 

of 

probability 


%Coordinates of emitter 
E=[0 0 0 ] ; 

%Speed of Light(Propagation) 
c = 3 * 10 A 8 ; 

%Set All Alphas to 1 

alphal = 1; 

alpha2 = 1; 

alpha3 = 1; 

alpha4 = 1; 

alpha5 = 1; 
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%Stationary Receivers 


X2 

= 1000; 

Y2 = 

-2000; 

Z2 = 

-2100; 

X3 

=-1000; 

Y3 = 

2000; 

Z3 = 

-2000; 

X4 

= -1000; 

Y4 = 

-1000; 

Z4 = 

1000; 

X5 

= -3000; 

Y5 = 

-1000; 

Z5 = 

-1000; 


%for loop for UAV movement 
for j = l:m, 

%locations of the UAV 

XI = 10 0 0 * j; Y1 = 15 0 0 * j; Z1 = 2000; %Moving UAV 


%UAV Coordinate MAtrix 
XlCoor(1,j) = XI; 

YlCoor (1,j) = Yl; 

ZlCoor(1,j) = Zl; 

%for loop for experiments for every location of the UAV 
for i = 1;n, 

%initial coordinate matrix 
Coor2 = [0 0 0]; 

%Calculate the Distance From emitter to the 
Receiver for TOA 

ECl=sqrt(abs( (E (1)-XI) A 2+(E(2)-Yl) A 2+(E(3)-Z1) A 2)) ; 
EC2=sqrt(abs((E(1)-X2) A 2+(E(2)-Y2) A 2+(E(3)-Z2) A 2)); 
EC3=sqrt(abs((E(1)-X3) A 2+(E(2)-Y3) A 2+(E(3)-Z3) A 2)); 
EC4=sqrt(abs( (E(1)-X4) A 2+(E(2)-Y4) A 2+(E(3)-Z4) A 2)) ; 
EC5=sqrt(abs((E(1)-X5) A 2+(E(2)-Y5) A 2+(E(3)-Z5) A 2)); 

%Generate random noise for every channel related to 
Standat Deviation 

Noise Err = random(' Normal 0,5e-8,1,5); 


%Calculate the times of arrival. 


tl 

= EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err(1) 

r 

t2 

= EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err (2; 

r 

t3 

= EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err(3] 

r 

t4 

= EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err(4) 

r 

t5 

= EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err (5) 

r 


%from now on calculate the emitter location 
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%Calculate the coefficients in Eqs.(17). Note that 
alpha (assumed),not alpha actual is used. 

all = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t3-tl) *(X3/alpha3 A 2 - XI/alphal A 2); 

al2 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t3-tl) * (Y3/alpha3 A 2 - Y1/alphal A 2); 

al3 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t3-tl) * (Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 

((X3 A 2+Y3 A 2+Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (t3-t2); 

a21 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/ (14-t1) *(X4/alpha4 A 2 - XI/alphal A 2); 

a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/ (14-t1) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t4-tl) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (14-12) ; 

a31 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t5-tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t5-tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t5-tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

((X5 A 2+Y5 A 2+Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (15-12) ; 

%Calculate the matrix A on the left-hand side. 

A = double([[all al2 al3]; [a21 a22 a23]; [a31 a32 

a 3 3 ] ] ) ; 


%Calculate the vector B on the right-hand side. 
B = double([bl; b2; b3]); 


%Solve the simultaneous equations AX = B. 
Coorl = (A'*A)\A'*B; 

%Error by every Coordinate 
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Err = [abs (Coorl (1)-E (1)),abs (Coorl (2)-E (2)), 
abs (Coorl (3)-E (3))]; 

%Error Matrix 
Mat(i,:) = Err; 


end 


%calculate Cov.Matrix 
Cov_Mat = cov(Mat); 

%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues ] = eig (Cov_M a t) ; 

%Find Diagonalized Cov.Mat. 

Diagonal_C° v _Mat = 

Eigen_V ec tors.'*inv(Cov_M & t)*Eigen_V ec tors; 

%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*k); 

Std_Dev_Data_Mat(:,j) = Std_Dev_Data 

%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 


Tot_Std_Dev_M a t (j,:) = Tot_Std_Dev 


end 

figure 

%Draw the error ellipsoid from the diagonalized covariance 

matrix 

hold on 

grid on 

%equal axis 
axis equal 

%# Change the camera viewpoint 
view([-36 18]); 

%label axises and title 
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xlabel( 'X-Axis' ) 
ylabel( 'Y-Axis' ) 
zlabel( 'Z-Axis' ) 
title( 'TDOA Geometry') 

%plot stationary receivers 

plot3 (X2,Y2,Z2, ' o ' ) ; 
plot3 (X3,Y3,Z3, 'o' ); 
plot3 (X4,Y4,Z4, 'o' ) ; 
plot3 (X5,Y5,Z5, 'o' ); 

%plot emitter 

plot3 (E (1) ,E (2) ,E (3) , 'ms ' ) ; 

%plot moving receiver 

plot3(XlCoor,YlCoor,ZlCoor, 'r—' ); 

figure 

N = linspace(1,m,m); 

%plot the result 

plot (N,Tot_Std_Dev_Mat) 

%label axises and title 

xlabel (' Distance Multiplier') 

ylabel (' Total Standard Deviation (m) ') 

title (' Effect of Distance to the Accuracy With UAV') 
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MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 
TECHNIQUE TO SEE THE EFFECTS OF HAVING TWO MOVING UAVS 
ALONG WITH STAIONARY RECEIVERS ON ACCURACY. 


function TDOA_Flying_two_UAV 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 

% Equations, 2006 

% Latest revision: August 26, 2012 

%clear the screen 
clc 

%Variables 

%number of experiments for accuracy 
n = 500; 

%number of the positions of the UAV 
m = 2 0; 


k=8.934259824; %k for 95% of probability 


o 

o 

k=6.922544695; 

%k 

for 

90% 

of 

probability 

o 

o 

k=5.766298681; 

%k 

for 

85% 

of 

probability 

o, 

o 

k=4.949459443; 

%k 

for 

80% 

of 

probability 

o, 

o 

k=4.314831079; 

%k 

for 

75% 

of 

probability 


%Coordinates of emitter 
E=[0 0 0 ] ; 

%Speed of Light(Propagation) 
c = 3 * 10 A 8 ; 


%Set All Alphas to 1 


alphal = 1; 
alpha2 = 1; 
alpha3 = 1; 
alpha4 = 1; 
alpha5 = 1; 




%Stationary 

Receivers 



X3 =-1000; 

Y3 = 2000; 

Z3 = 

-2000; 

X4 = -1000; 

Y4 = -1000; 

Z4 = 

1000; 

X5 = -3000; 

Y5 = -1000; 

Z5 = 

-1000; 
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%for loop for UAV movement 
for j = l:m, 

%locations of the UAV 

XI = 8 0 0 * j; Y1 = 12 0 0 * j; Z1 = 2000; %Moving UAV 

X2 = 7 0 0 * j; Y2 = -1800*j; Z2 = 2000; 

%UAV Coordinate MAtrix 
XlCoor(1,j) = XI; 

YlCoor (1,j) = Yl; 

ZlCoor(1,j) = Zl; 

X2Coor (1,j) = X2; 

Y2Coor(1,j) = Y2; 

Z2Coor(1,j) = Z2; 

%for loop for experiments for each location of the UAV 
for i = 1:n, 

%Calculate the Distance From emitter to the 
Receiver for TOA 

ECl=sqrt(abs( (E(1) -XI) A 2+(E(2)-Yl) A 2+ (E (3)-Zl) A 2)) ; 
EC2=sqrt(abs( (E(1)-X2) A 2+(E (2)-Y2) A 2+ (E (3)-Z2) A 2)) ; 
EC3=sqrt(abs((E(1)-X3) A 2+(E(2)-Y3) A 2+(E(3)-Z3) A 2)); 
EC4=sqrt (abs ( (E (1) -X4) A 2+ (E (2) -Y4) A 2+(E(3)-Z4) A 2)); 
EC5=sqrt (abs ( (E (1) -X5) A 2+ (E (2) -Y5) A 2+(E(3)-Z5) A 2)); 

%Generate random noise for every channel related to 
Standat Deviation 

Noise Err = random(' Normal 0,5e-8,1,5); 


%Calculate 

the times 

of 

arrival 

. 


tl 

= EC1 

/ 

(alphal 

* 

c) 

+ 

Noise 

Err 

(l) ; 

t2 

= EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err(2) 

t3 

= EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

Err(3 ); 

14 

= EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err 

(4) ; 

t5 

= EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err 

(5) ; 


%from now on calculate the emitter location 

%Calculate the coefficients in Eqs.(17). Note that 
alpha (assumed),not alpha actual is used. 

all = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t3-tl) *(X3/alpha3 A 2 - XI/alphal A 2); 
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al2 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t3-tl) *(Y3/alpha3 A 2 - Y1/alphal A 2); 

al3 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 

2/(t3-tl) *(Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 

((X3 A 2+Y3 A 2+Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (13-12); 

a21 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/ (14-t1) *(X4/alpha4 A 2 - XI/alphal A 2); 

a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t4-tl) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/ (14-t1) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (14-12) ; 

a31 = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 
2/(t5-tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 
2/(t5-tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 
2/(t5-tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/(t2-tl) * ((X2 A 2+Y2 A 2+Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

( (X5 A 2+Y5 A 2 + Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2 + Z1 A 2)/alphal A 2) + 
c A 2 * (15-12) ; 

%Calculate the matrix A on the left-hand side. 

A = double([[all al2 al3]; [a21 a22 a23]; [a31 a32 

a 3 3 ] ] ) ; 


%Calculate the vector B on the right-hand side. 
B = double([bl; b2; b3]); 

%Solve the simultaneous equations AX = B. 

Coorl = (A'*A)\A'*B; 

%Error by every Coordinate 

Err = [abs(Coorl(1)-E (1)),abs (Coorl (2)-E (2)), 
abs (Coorl (3)-E (3))]; 

%Error Matrix 
Mat (i, :) = Err; 


134 


end 


%calculate Cov.Matrix 
Cov_Mat = cov(Mat); 

%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues] = eig (Cov_M a t) ; 

%Find Diagonalized Cov.Mat. 

Diagonal_C° v _^ a t = 

Eigenvectors _ 1 *j_nv (Cov_M a t) *Eigen_V ec tors; 

%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*k); 

Std_Dev_Data_M a t(:,j) = Std_Dev_Data; 

%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 


Tot_Std_Dev_M a t (j,:) = Tot_Std_Dev 


end 

figure 

%Draw the error ellipsoid from the diagonalized covariance 

matrix 

hold on 

grid on 

%equal axis 
axis equal 

%# Change the camera viewpoint 
view([-36 18]); 

%label axises and title 
xlabel( 'X-Axis' ) 
ylabel( 'Y-Axis' ) 
zlabel( 'Z-Axis' ) 
title( 'TDOA Geometry') 
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%plot stationary receivers 


plot3 (X3,Y3,Z3, 'o' ); 
plot3 (X4,Y4,Z4, 'o' ) ; 
plot3 (X5,Y5,Z5, 'o' ); 

%plot emitter 

plot3 (E (1) ,E (2) ,E (3) , 'ms ' ) ; 

%plot moving receiver 

plot3(XlCoor,YlCoor,ZlCoor, 'r--' ); 

plot3(X2Coor,Y2Coor,Z2Coor, 'g—' ); 

figure 

N = linspace(1,m,m); 

%plot the result 

plot (N,Tot_Std_Dev_Mat) 

%label axises and title 

xlabel (' Distance Multiplier') 

ylabel (' Total Standard Deviation (m) ') 

title (' Effect of Distance to the Accuracy With UAV') 
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o\° o\° 


APPENDIX G 


MATLAB CODE FOR SIMULATION OF CLOSED-FORM GEOLOCATION 
TECHNIQUE TO SEE THE EFFECTS OF THE RECEIVERS GEOMETRY ON 

ACCURACY. 

function TDOA_Receiver_Dist 
% Written by: Volkan TAS 

% Based on code developed by EZZAT G. BAKHOUM, described in 
% Closed-Form Solution of Hyperbolic Geolocation 
% Equations, % 2006 
% Latest revision: August 15, 2012 

%clear the screen 
clc 

%Variables 

%number of experiments for accuracy 
n = 500; 

Mat = zeros(n,3); 


k=8.934259824; %k for 95% of probability 


% k=6.922544695; 

%k 

for 

90% 

of 

probability 

% k=5.766298681; 

%k 

for 

85% 

of 

probability 

% k=4.949459443; 

%k 

for 

80% 

of 

probability 

% k=4.314831079; 

%k 

for 

75% 

of 

probability 

%Coordinates of 

emitter 




E=[0 0 0 ] ; 






%Speed of Light(Propagation) 




c = 3 * 10 A 8 ; 

%Set All Alphas to 1 

alphal = 1; 

alpha2 = 1; 

alpha3 = 1; 

alpha4 = 1; 

alpha5 = 1; 

Stationary Receivers 
set default for switch 
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val 


1 ; 


tmp = input (sprintf ('How you want to Distribute the 
Receivers \n Around the Emitter (Without Altitude 
Difference)Press 1 \n On Staright Line (Without Altitude 
Difference)Press 2 \n Around the Emitter (With Altitude 
Difference) Press 3 \n On Staright Line (With Altitude 
Difference)Press 4 \n Trapezoidal Press 5 \n 
Parallelogram Press 6 \n Inverted Triangle Press 8 
' , 'val' )); 

%# if the user hits 'return' without writing anything, tmp 
is empty and the default is used 
if -isempty(tmp) 
val = tmp; 

end 

switch val 


case 1 

%circle, no 

elevation difference 


XI 

= 2 e 3 ; 

Y1 = 

0; 

Zl = 

1.9 7 e 3 ; 

X2 

= 7e2 ; 

Y2 = 

2e3; 

Z2 = 

1.9 8 e 3 ; 

X3 

=-2e3; 

Y3 = 

1.5 e 3; 

Z3 = 

1.96e3; 

X4 

= -2e3; 

Y4 = 

-1.5 e 3; 

Z4 = 

2e3; 

X5 

= 1 e 3 ; 

Y5 = 

-1.5e3; 

Z5 = 

1.99e3; 

case 2 

%0n Staright 

Line 

(Without 

Altitude Difference 

XI 

= -2e3; 

Y1 = 

-1 e 3 ; 

Zl = 

0.99e3; 

X2 

= -1 e 3 ; 

Y2 = 

0; 

Z2 = 

le3 ; 

X3 

= 0; 

Y3 = 

2e3 ; 

Z3 = 

0.9 8 e 3 ; 

X4 

= 1 e 3 ; 

Y4 = 

3e3 ; 

Z4 = 

0.97e3; 

X5 

= 1.5 e 3; 

Y5 = 

4e3; 

Z5 = 

1.Ie3; 

case 3 

%Around the 

Emitter (With 

Altitude Difference) 

XI 

= 2 e 3 ; 

Y1 = 

0; 

Zl = 

0; 

X2 

= 0.7 e 3 ; 

Y2 = 

2e3 ; 

Z2 = 

le3 ; 

X3 

=-2e3; 

Y3 = 

1.5e3; 

Z3 = 

-1 e 3 ; 

X4 

= -2e3; 

Y4 = 

-1.5e3; 

Z4 = 

2e3 ; 

X5 

= le3 ; 

Y5 = 

-1.5 e 3; 

Z5 = 

-2e3 ; 

case 4 

%On Staright 

Line 

(With Altitude 

Difference) 

XI 

= -2e3; 

Y1 = 

-1 e 3 ; 

Zl = 

0; 

X2 

= -1 e 3 ; 

Y2 = 

0; 

Z2 = 

le3 ; 

X3 

= 0; 

Y3 = 

2e3 ; 

Z3 = 

2e3 ; 

X4 

= le3 ; 

Y4 = 

3e3 ; 

Z4 = 

3e3 ; 

X5 

= 1.5e3; 

Y5 = 

4e3 ; 

Z5 = 

4e3; 
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case 5 %Trapezoidal 


XI 

= le3 ; 

Yl 

= 

0; 

Zl 

= 

0; 

X2 

= 2 e 3 ; 

Y2 

= 

2e3 ; 

Z2 

= 

-0.1e3; 

X3 

= 3e3; 

Y3 

= 

2e3 ; 

Z3 

= 

-0.1e3; 

X4 

= 4e3 ; 

Y4 

= 

2e3 ; 

Z4 

= 

- 0.1 e 3 ; 

X5 

= 5e3; 

Y5 

= 

0; 

Z5 

= 

-0.1e3; 

case 6 

%Parallelogram 






XI 

= -1 e 3 ; 

Yl 

= 

0; 

Zl 

= 

0; 

X2 

= -2e3; 

Y2 

= 

2e3 ; 

Z2 

= 

- 0.1 e 3; 

X3 

= 0; 

Y3 

= 

2e3 ; 

Z3 

= 

- 0.1 e 3; 

X4 

= -2e3; 

Y4 

= 

0; 

Z4 

= 

- 0.1 e 3; 

X5 

= -3e3; 

Y5 

= 

0; 

Z5 

= 

- 0.1 e 3 ; 

case 7 

%Lorenze 







XI 

= 8e3; 

Yl 

= 

15e3; 

Zl 

= 

- 0.12 e 3 ; 

X2 

= -13 e 3 ; 

Y2 

= 

12.5e3; 

Z2 

= 

-0.13e3; 

X3 

= 0; 

Y3 

= 

18 e 3 ; 

Z3 

= 

-0.14e3; 

X4 

= 13e3; 

Y4 

= 

12.5 e 3; 

Z4 

= 

- 0.11 e 3; 

X5 

= 0; 

Y5 

= 

8; 

Z5 

= 

- 0.1 e 3; 

case 8 

%Inverted 

Triangle 




XI 

= 2 e 3 ; 

Yl 

= 

0; 

Zl 

= 

0; 

X2 

= 4 e 3 ; 

Y2 

= 

2e3 ; 

Z2 

= 

- 0.1 e 3; 

X3 

= 6e3; 

Y3 

= 

-2e3; 

Z3 

= 

-0.1e3; 

X4 

= 8e3; 

Y4 

= 

4e3; 

Z4 

= 

-0.1e3; 

X5 

= 10 e 3; 

Y5 

= 

- 4 e 3; 

Z5 

= 

- 0.1 e 3; 


end 

for i = 1:n, 

%Calculate the Distance From emitter to the Receiver 
for TOA 

ECl=sqrt(abs( (E(1)-XI) A 2+(E(2)-Yl) A 2 +(E (3)-Zl) A 2)) ; 
EC2=sqrt(abs( (E(1)-X2) A 2+(E(2)-Y2) A 2+ (E (3)-Z2) A 2)) ; 
EC3=sqrt(abs((E(1)-X3) A 2+(E(2)-Y3) A 2+(E(3)-Z3) A 2)); 
EC4=sqrt(abs( (E(1)-X4) A 2+ (E (2)-Y4) A 2+ (E (3)-Z4) A 2)) ; 
EC5=sqrt(abs( (E(1)-X5) A 2+(E(2)-Y5) A 2+ (E (3)-Z5) A 2)) ; 

%Generate random noise for every channel related to 
Standat Deviation 

Noise_Err = random (' Normal 0,5e-8,1,5); 

%Calculate the times of arrival. 

tl = EC1 / (alphal * c) + Noise_Err(1); 
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t2 

= EC2 

/ 

(alpha2 

* 

c) 

+ 

Noise 

Err (2) ; 

t3 

= EC3 

/ 

(alpha3 

* 

c) 

+ 

Noise 

_Err (3) ; 

14 

= EC4 

/ 

(alpha4 

* 

c) 

+ 

Noise 

Err(4); 

t5 

= EC5 

/ 

(alpha5 

* 

c) 

+ 

Noise 

Err(5); 


%from now on calculate the emitter location 


%Calculate the coefficients in Eqs.(17). Note that 
alpha (assumed),not alpha actual is used. 

all = 2/(t2-tl) * (X2/alpha2 A 2 - Xl/alphal A 2) - 2/(t3 
tl) *(X3/alpha3 A 2 - XI/alphal A 2); 

al2 = 2/ (t2-11) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 2/(t3 
tl) *(Y3/alpha3 A 2 - Y1/alphal A 2); 

al3 = 2/ (t2-11) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 2/(t3 
tl) *(Z3/alpha3 A 2 - Zl/alphal A 2); 

bl = 1/(t2-tl) * ( (X2 A 2+Y2 A 2 + Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t3-tl) * 

( (X3 A 2+Y3 A 2 + Z3 A 2)/alpha3 A 2- (X1 A 2+Y1 A 2 + Z1 A 2)/alphal A 2) + 
c A 2 * (t3-t2) ; 

a21 = 2/ (t2-11) * (X2/alpha2 A 2 - Xl/alphal A 2) - 2/(t4 
tl) *(X4/alpha4 A 2 - XI/alphal A 2); 

a22 = 2/(t2-tl) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 2/(t4 
tl) *(Y4/alpha4 A 2 - Y1/alphal A 2); 

a23 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 2/(t4 
tl) *(Z4/alpha4 A 2 - Zl/alphal A 2); 

b2 = 1/ (t2-t1) * ((X2 A 2+Y2 A 2 + Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t4-tl) * 

((X4 A 2+Y4 A 2+Z4 A 2)/alpha4 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (14-12) ; 

a31 = 2/ (t2-11) * (X2/alpha2 A 2 - Xl/alphal A 2) - 2/(t5 
tl) *(X5/alpha5 A 2 - XI/alphal A 2); 

a32 = 2/ (t2-11) * (Y2/alpha2 A 2 - Yl/alphal A 2) - 2/(t5 
tl) *(Y5/alpha5 A 2 - Y1/alphal A 2); 

a33 = 2/(t2-tl) * (Z2/alpha2 A 2 - Zl/alphal A 2) - 2/(t5 
tl) *(Z5/alpha5 A 2 - Zl/alphal A 2); 

b3 = 1/ (t2-t1) * ((X2 A 2+Y2 A 2 + Z2 A 2)/alpha2 A 2- 

(XI A 2+Y1 A 2+Z1 A 2)/alphal A 2)- l/(t5-tl) * 

((X5 A 2+Y5 A 2+Z5 A 2)/alpha5 A 2- (X1 A 2+Y1 A 2+Z1 A 2)/alphal A 2) + 
c A 2 * (t5-t2) ; 

%Calculate the matrix A on the left-hand side. 

A = double([[all al2 al3]; [a21 a22 a23]; [a31 a32 

a 3 3 ] ] ) ; 
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%Calculate the vector B on the right-hand side. 
B = double([bl; b2; b3]); 

%Solve the simultaneous equations AX = B. 

Coor = (A'*A)\A'*B; 


%Error by every Coordinate 

Err = [abs(Coor (1)-E(1)),abs (Coor (2)-E (2)), 
abs (Coor (3)-E (3))]; 

%Error Matrix 
Mat (i, :) = Err; 


end 

%calculate Cov.Matrix 
Cov_Mat = cov(Mat); 

%Eigenvalues and eigenvectors 
[Eigenvectors, Eigenvalues ] = eig(Cov_M a t); 

%Find angles (from "Computing Euler angles from a rotation 
matrix") 

Beta = -asind (Eigenvectors (3, 1) ) ; 

Gamma = 

(atan2 (Eigen_V ec tors(2,1)/cosd(Beta),Eigen_V ec tors(1, 1) / cos 
d(Beta)))*(180/pi); 

Alpha = 

(atan2(Eigen_V ec tors(3,2)/cosd(Beta),Eigen_V ec tors(3,3)/cos 
d(Beta)))*(180/pi); 

%Std. Dev (includes k value) 

Std_Dev_Data = sqrt(diag(Eigenvalues)*k) 

%Total Standat Deviation 
Tot_Std_Dev = 

sqrt(Std_Dev_Data(1) A 2+Std_Dev_Data(2) A 2+Std_Dev_Data(3) A 2) 

%Draw the error ellipsoid from the diagonalized covariance 

matrix 

hold on 

grid on 
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%ellipsoid 
[x, y, z] = 

ellipsoid(E(1),E(2),E(3),Std_Dev_Data(1),Std_Dev_Data(2),St 
d_Dev_Data(3),20); 

%color of the ellisoid 
colormap copper; 

%ploy ellipsoid 
hMesh = mesh(x,y,z); 

%rotate ellipsoid 
rotate(hMesh,[1 0 0],Alpha); 
rotate(hMesh,[0 1 0],Beta); 
rotate(hMesh,[0 0 1],Gamma); 

%equal axis 
axis equal 

%# Change the camera viewpoint 
view([-36 18]); 

%label axises and title 
xlabel( 'X-Axis' ) 
ylabel( 'Y-Axis' ) 
zlabel( 'Z-Axis' ) 

title( 'TDOA Geometry and Confidence Ellipsoid') 

%plot stationary receivers 
plot3 (XI,Yl,Zl, 'o' ) ; 
plot3 (X2,Y2,Z2, 'o' ); 
plot3 (X3,Y3,Z3, 'o' ); 
plot3 (X4,Y4,Z4, 'o' ); 
plot3 (X5,Y5,Z5, 'o' ); 

%plot emitter 

plot3 (E (1),E (2),E (3), 'ms' ); 
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